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A coupled Eulerian/Lagrangian method is presented for the reduction of numerical 
diffusion observed in solutions of three-dimensional rotational flows using standard Eu- 
lerian finite-volume time-marching procedures. A Lagrangian particle tracking method 
using particle markers is added to the Eulerian time-marching procedure and provides 
a correction of the Eulerian solution. In turn, the Eulerian solution is used to integrate 
the Lagrangian state- vector along the particles trajectories. The Lagrangian correction 
technique does not require any a-priori information on the structure or position of the 
vortical regions. While the Eulerian solution ensures the conservation of mass and sets 
the pressure field, the particle markers, used as ‘accuracy boosters’, take advantage of 
the accurate convection description of the Lagrangian solution and enhance the vorticity 
and entropy capturing capabilities of standard Eulerian finite-volume methods. 


The combined solution procedure is tested in several applications. The convection of 
a Lamb vortex in a straight channel is used as an unsteady compressible flow preservation 
test case. The other test cases concern steady incompressible flow calculations and 
include the preservation of a turbulent inlet velocity profile, the swirling flow in a pipe, 
the constant stagnation pressure flow and secondary flow calculations in bends. The 
last application deals with the external flow past a wing with emphasis on the trailing 
vortex solution. 


The improvement due to the addition of the Lagrangian correction technique is mea- 
sured by comparison with analytical solutions when available or with Eulerian solutions 
on finer grids. The use of the combined Eulerian/Lagrangian scheme results in substan- 
tially lower grid resolution requirements than the standard Eulerian scheme for a given 
solution accuracy. 
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Chapter 1 


Introduction 


1.1 Statement of the problem 


Over the last few years, the improvement in CPU and memory capabilities of modem 
supercomputers has rendered practical the solution of flow problems of more and more 
complex nature. However, the efficient numerical treatment of flow non-homogeneities, 
such as vortex wakes or tip vortex roll-up, embedded in an otherwise smooth background 
flow field remains a challenging field of study. In many practical applications, the 
prediction of the strength and the position of the vortical regions reveals to be of primary 
importance. For instance, the flow around am helicopter rotor blade presents a case of 
strong interaction between the shed vortices due to one blade and the following blade. 
The prediction of the resulting load variations requires the accurate solution of the shed 
vortices, of their trajectories and of the subsequent interaction phenomenon. Another 
example is the prediction of the secondary flow through a bend with a pump attached 
at the bend exit. The location and strength of the secondary vortex, created by the 
tilting and stretching of the inlet boundary-layer vorticity, must be solved accurately, 
as a possible noise source and a performance loss may result from the impingement of 
the secondary vortex on the rotating blades. 

Vortex- sheets, secondary flows, or vortex roll-up phenomena are all characterized 
by transverse length scales differing by orders of magnitude from the length scale of 
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the supporting flow field (the transverse length scale of trailing vortices has been found 
experimentally to be as low as 5% of the airfoil chord [57]). Since the solution of 
these vortical features has often to cover a convection length much higher than their 
intrinsic length scale, the global prediction of vortex-dominated flows proves to be highly 
sensitive to sm all local errors. This makes these flow features difficult to be captured 
by finite-difference methods. 


1.2 Existing approaches 


Incompressible vortex methods and potential methods with fitted vortex sheets are 
not susceptible to the numerical diffusion. Examples of incompressible vortex methods, 
using Biot-Savart law to compute the velocity field, include the method used by Leonard 
[41] where the flow vorticity is modeled as a collection of a few isolated vortex-tubes 
with a computational element assigned to each vortex tube and Knio’s study [38] where 
a three-dimensional vortex scheme is based on the transport of vorticity and material 
elements. Additionally, gradients of the scalar field are transported and the scalar field 
itself is recovered using Biot-Savart law. 

Potential methods presuppose some a priori knowledge (either from a known solution 
or from empirical data) on the vortex structure or position, since potential methods 
do not ‘capture’ embedded vorticity as part of the solution. This limitation becomes 
especially acute when solving complex flow problems where an a priori information is 
not always available. Scully [64] and Miller [44] have used a Biot-Savart formulation for 
an incompressible flow solution of an helicopter rotor wake. Hassan [25] used the Euler 
equations in an implicit scheme and modeled the blade- vortex interaction by computing 
the vortex-induced velocities following Biot-Savart law. Steinhoff [71, 72] presented an 
alternative method for an aircraft configuration, where the strength, position and shape 
of the vortex sheet were calculated as part of the solution. The internal structure of the 
wake was, however, still to be specified. Ramach an dr an [54] used a potential method 
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with embedded vortex wakes for the compressible flow solution of a rotor wake. The 
body was included in the calculation similarly as the rotor wake as a vorticity sheet 
whose strength is determined iteratively. 

The use of an Euler or a Navier- Stokes solver presents the advantage that the em- 
bedded vorticity is captured as part of the solution. Poor solution representation is, 
however, a common feature of Eulerian and Navier- Stokes solvers in regions of high gra- 
dients. More precisely, the errors introduced by discretizing the equations of motion can 
be expressed in terms of dissipation and dispersion phenomena. In addition, rounding 
errors are introduced randomly in the solution. The dissipation expresses the fact that 
the finite difference model loses energy as the time progresses. Because numerical dis- 
sipation can be advantageous by counteracting unwanted instabilities and oscillations 
(for example saw-tooth modes), it is added to non- dissipative formulae. The disper- 
sion errors correspond to the decay of a wave form into separate spurious oscillations 
and always occur with finite-difference formulae since their dispersion relation is always 
non-linear. 

Because most common second-order accurate finite-difference schemes will smear 
and distort regions of high gradients, corresponding grid clustering seems the obvious 
approach. However, depending on the flow topology, this procedure could reveal to 
be prohibitively expensive. For example, the solution of the flow around helicopter 
blades involves the prediction of the interaction between the tip vortex from a blade 
and the following blade. Because the resulting flow presents vortex regions of high 
extent and is highly non-homogeneous, standard clustering of the grid would lead to 
high computational cost. As Drela [20] reported, a suitably overall fine grid would imply 
~ 260 billion points for a rotor wake solution. The sensitivity of the solution to the grid 
coarseness has, therefore, prompted several studies. 

Selective refinement of the grid is used by Lohner [42] in an adaptive algorithm 
in two dimensions. In three dimensions, however, the algorithm would undoubtedly 
present much complexity. Nakahashi [46] uses a 2-D solution adaptive structured grid 


21 



method based on variational principles and spring analogies. A multigrid solution of 
the Euler equations using an implicit scheme has been performed by Jameson [37] and 
has proven effective for CPU reduction in two-dimensional applications. The three- 
dimensional application would certainly prove to be more involved. Landsberg [40] 
has studied vortex capturing using adaptation and a three-dimensional finite-element 
solver. Also Schmatz [63] presented a two-dimensional zonal solution to model the weak 
or strong viscous/inviscid interactions in subsonic and transonic flows. Powell [49] used 
an adaptive mesh procedure working on an unstructured mesh for solving the conical 
Euler equations for leading-edge vortex flows. In this respect, unstructured grids present 
a strong advantage over structured grids because of the flexibility of adding new grid 
resolution in defined areas. 

An alternative to grid refinement is to use a high-order accurate scheme, a method 
used by Rai [53] who presented a fifth-order upwind-biased scheme in a blade/ vortex 
interaction problem. Steger [70] used an implicit fourth-order accurate scheme for the 
computation of vortex wakes. However, the advantage of using these high-order accurate 
schemes, with suitable clustering of points in the regions of high gradients, is still linked 
to grid smoothness (more difficult to obtain in three dimensions) and were demonstrated 
on grids prohibitively fine for complex three-dimensional applications. 

Perturbation methods, instead, rely on a known 1 flow solution in some areas of the 
computational domain such as to correct the numerical diffusion encountered in the 
basic finite- difference flow solution. Roberts [57, 56] applied this methodology, first 
introduced by Chow [15] to the rotor wake and blade/ vortex interaction problems by 
coupling the Euler equations to a free-wake model of the rotary wing wake. The main 
drawback here is the need for a known solution which limits the correction method 
to regions of simple behavior (the correction can not be performed where the vortex 
impinges on the airfoil for instance). In a similar approach, Srinivasan [69, 68] uses the 
2-D thin-layer Navier-Stokes equations and a prescribed vortex for computing the flow 

1 This is done from an analytic solution if available or from a previously computed high resolution 
local solution of the flow field. 
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over rotor blades. 


The ‘Cloud in Cell’ technique, first introduced by Christiansen [16] and then trans- 
formed by Baker [5] for a 2-D incompressible inviscid fluid, uses an area averaged vor- 
ticity distribution from markers in cells onto grid nodes where Poisson’s equation is 
solved. The circulation distribution must, however, be known at each marker and a 
three- dimensioned solution is not straightforward. Also Basuki [7] used the inviscid 
‘Cloud in Cell’ technique with vortices tracked though the grid on which the velocity is 
found by a finite- difference method. Poor resolution of the velocity field was, however, 
reported. 

The advantages of spectral methods are accuracy, ease of implementation and the 
low number of collocation points required for a computation when compared to the dis- 
cretization used in finite- difference methods. However, they lead to large matrices for 
more than one non-periodical direction (which is the case for the flow cases treated here) 
and are difficult to apply in computational domains of complex shapes. Also, no discon- 
tinuity (as a shock, for example) is allowed as part of the solution. Furthermore they 
are more time restrictive than finite- difference methods for unsteady flow calculations. 


1.3 Present approach 


As mentioned in the previous section, the Euler or Navier- Stokes equations present 
the advantage of directly capturing embedded vorticity, when compared to potential 
methods with fitted sheets. For example, as reported by Murman [45] and others [50, 9], 
the use of the Euler equations provides a good tool for a study of flows around wings 
enabling the study of the leading-edge vortex. The capturing capabilities of Euler 
or Navier-Stokes solvers are, nevertheless, limited by grid resolution issues. As an 
alternative to the methods dealing with this problem, this thesis presents a technique 
for substantially improving the capturing capabilities of time-marching Eulerian solvers. 
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This method does not require the knowledge of an a priori solution, grid refinement or 
the use of higher- order schemes. 

The objective of this research is then to construct an alternative solution procedure to 
reduce the numerical diffusion observed in standard Eulerian time-marching calculations 
and to demonstrate the feasibility, efficiency and flexibility of the method by application 
to different flow problems. These include three-dimensional steady, unsteady, internal as 
well as compressible and incompressible inviscid flow cases. Furthermore, the extension 
of the method to include the Navier-Stokes equations (not performed in the frame of 
this thesis) is judged to be straightforward. 

The present method consists in the addition of a Lagrangian particle tracking so- 
lution to a standard Eulerian solution in order to enhance its vorticity and entropy 
capturing capabilities. This method is based on the approach of Drela [20] in two di- 
mensions and is here extended to include three-dimensional flow cases. The combination 
of the Eulerian and Lagrangian solvers takes advantage of both the accurate convection 
description of the Lagrangian technique and the ‘elliptic’ representation of the Eulerian 
solution which enforces the mass conservation and sets the pressure field. Briefly, the 
Lagrangian solution is based on particle markers carrying vorticity and entropy, and 
convecting with the local flow through the Eulerian grid. The Eulerian solver is used 
to conserve mass and to provide the source terms required for the Lagrangian time 
integration. In turn, since the Lagrangian solution is immune to the numerical diffu- 
sion process occurring in the Eulerian solver, it accurately captures the convection of 
vorticity and entropy. This information is then used to locally correct the Eulerian so- 
lution and to reduce its numerical diffusion errors. Each Lagrangian marker influences 
the Eulerian solution only locally (as opposed to the ‘Cloud in Cell 5 technique where 
each marker has an influence on the entire flow field), which makes this scheme well 
suited for three-dimensional flow solutions. Also, the Lagrangian solution needs to be 
computed only in regions of interest as markers can be located selectively in the flow. 
No a priori information is required on the flow structure since the Lagrangian solution 
includes inherently ‘convective 5 capabilities. 
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1.4 Thesis outline 


Chapters 2 and 3 present the equations, numerical procedure and accuracy study for 
the Eulerian solver for both the compressible and incompressible flow cases. The mesh 
generation technique is described in Appendix A. The Lagrangian equations are the 
object of Chapter 4, whereas the coupling of the two solvers in the time integration is 
described in Chapter 5 for different flow configurations. Finally, the correction proce- 
dure by which the Lagrangian equations influence the Eulerian solution is presented in 
Chapter 6. 

The first test case, discussed in Chapter 7, is the convection of a Lamb vortex in 
a three-dimensional uniform background flow and is used as a preservation test for a 
compressible unsteady flow. An analogous test case is presented in Chapter 8 as the 
preservation of a turbulent inlet velocity profile in a straight pipe. The swirling flow in 
a straight pipe is the object of Chapter 9 with an emphasis on the development of a 
vorticity gradient augmentation phenomenon and the particular solution adopted with 
the combined Eulerian /Lagrangian solver. The vorticity errors and stagnation pressure 
losses encountered in the Eulerian solution of a 90° bend are reduced by the use of the 
Lagrangian correction method in Chapter 10. 

Chapter 11 deads with the secondary flow in bent pipes. The secondary flow genesis 
is first described and Eulerian and Eulerian/Lagrangian solutions are computed and 
compared with experiments. The introduction of a simple Taw of the wall’ model is 
attempted in order to deal with viscosity effects. 

The last case, reported in Chapter 12, is the external flow over a three-dimensional 
wing. This chapter emphasizes the spurious numerical diffusion of the tip vortex behind 
the trailing edge and the correction obtained using the combined Eulerian/Lagrangian 
scheme. Comparison with experimental data is also performed. 

Finally, Chapter 13 presents a summary, the contributions as well as the conclusions 
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and the recommendations for future work. 
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Chapter 2 


Eulerian Governing Equations 


Provided the tangential forces applied on fluid particles are small compared to the 
pressure forces, the fluid can be treated as inviscid. The evolution of an inviscid flow 
in time and space is described by the Euler equations [60, 3]. Here, these equations 
are used for the solution of both steady incompressible and unsteady compressible flow 
fields. The numerical solution procedure using a Lax-WendrofF scheme is the object of 
the next chapter. 


2.1 Euler equations for compressible flow 


For the solution of unsteady compressible flows, the solution is marched forward in 
time from an initial condition. The Euler equations expressed in a (sc,t/, z) Cartesian 
coordinates system and in conservation form are 

dU _ d£ 3G dE_ 

dt dx + dy + dz 


( 2 . 1 ) 
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2.2 Euler equations for incompressible flow 


For an incompressible unsteady flow, the state vector U and the fluxes F, G and H are 
written as 



where p is a constant and the ratio p/p is denoted p*. 
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The main problem in solving the Euler equations for an incompressible flow is to link 
the velocity changes to the pressure changes in a way that enforces the divergence-free 
condition. In two dimensions, the solution of steady or unsteady incompressible flows 
can be achieved through the stream function- vorticity formulation (therefore eliminating 
the pressure from the governing equations). In three-dimensional flow calculations, 
however, this solution technique becomes more complex and other solution procedures 
are usually sought. The vorticity- velocity formulation used by Dennis [17] replaces the 
two-dimensional stream function- vorticity formulation. The Poisson’s equation method 
developed by Harlow [24] consists of iteratively adjusting the pressure field by solving 
a Poisson type equation for the pressure change. Poisson’s equation is obtained from 
the requirement that the continuity equation must be satisfied. Since this method 
involves an iterative procedure it is, however, very time consuming for three-dimensional 
applications. 

A well-known class of solution procedures for steady compressible flows is the time- 
marching method where the full unsteady Euler equations are used and the solution 
evolves through a pseudo-unsteady process from an initial guess to the final steady- 
state. Nevertheless, in the limit of an incompressible flow, sound waves with very large 
speed tend to make the system stiff and render this inefficient. A well-known solution 
to this problem, and the method used in this work, is the artificial compressibility 
concept introduced by Chorin [14]. The purpose of this technique is to transform the 
character of the Euler equations for an incompressible flow from elliptic to hyperbolic 
by adding a time- dependent term in the continuity equation. This particular method 
has been successfully tested on an extensive set of internal and external incompressible 
flow problems [58, 13, 55, 61, 79]. 

To introduce a time- derivative of the pressure in the continuity equation, the di- 
vergence term is multiplied by the “artificial compressibility” parameter c\ so that the 
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modified state and flux vectors are defined as 


f jr \ 


u = 


U 

V 

w 


F = 


t clu 


U 2 +p * 

uv 

uw 


G = 


/ civ 


uv 


V 2 + p* 


vw 


H = 


l clw \ 


uw 


vw 


\ w 2 + p* / 


• (2.5) 


The value of the artificial compressibility parameter can be adjusted to increase the 
convergence rate of the time-marching procedure. When steady-state is reached, the 
modified system of equations reduces to the standard Euler equations for steady flow. 
Also pseudo-time stepping can be used since the unsteady process is of no interest here 
and only the steady-state solution is retained. 


The introduction of the artificial compressibility parameter results in giving finite 
speed to the propagating waves, in contrast to truly incompressible flow where the waves 
move with an infinite speed. The pseudo-speed of sound 0 is computed in Section 3.4 by 
analyzing the linearized Euler equations with 1-D variations. 0 depends on the artificial 
compressibility parameter c\ as 

0 = \Ju 2 + c\. ( 2 . 6 ) 


Chang [13] estimates the relation between the parameter c\ and the rate of convergence 
by loo king at the speed of the propagating waves. The time taken by a wave to travel 
from the inlet of the computational domain to the exit over a distance L and back is 


L L _ 2(3 L 
0 + u + 0-u c\ 


(2.7) 


This value represents the minimum time needed for convergence. If the time-step allowed 
for stability is A t, then N the number of time-steps required becomes 


N = 


2 0L 
cl&t 



( 2 . 8 ) 


a decreasing function of c\. Regarding the value of the artificial compressibility param- 
eter, it is shown in Section 3.4 that the ratio of the largest to the smallest eigenvalues 
of the linearized Euler system of equations is dependent on 0. 0 is therefore a measure 


30 



of the condition of the system. Rizzi [55] has verified numerically that a ratio c 2 /u 2 
between 1 and 5 ensures the system to be well- conditioned. In this particular study a 
constant value of c\ju 2 of 1 is used. 

Another advantage of the artificial compressibility method is its natural extension 
to the handling of the Navier-Stokes equations [13]. Also the same concept has been 
used for the solution of unsteady flows problems [43]. 


2.3 Non-dimensionalization 


The Euler equations are used in a non-dimensionalized form which allows for the flow 
values to fall within prescribed limits. The arbitrary reference values are given below 
for the different flow cases treated in this work. 


For the compressible flow in a channel the reference quantities are the channel length 
Z, the inlet stagnation speed of sound cq iu and the inlet stagnation density po in so that 
the non-dimensional variables are 

= X’ b Zr= X’ 
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The corresponding reference inlet stagnation pressure and enthalpy are 


iPO tn )r • — 



1 

7 


(ho in )r 


7 (PO <n )r 
7-1 iPOiJr 



(2.9) 


For incompressible flow in pipes, the reference quantities are the pipe radius R, and 
the inlet mass-flow averaged velocity tZ^\ For the incompressible flow around an airfoil 
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the reference length has been chosen as the airfoil chord Ch and the reference velocity 
is the freestream velocity Uoo- If the reference length and the reference velocity are 
denoted by L re f and U Te j, then 
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Chapter 3 


Euler Solver: Numerical Procedure 


3.1 Lax-Wendroff algorithm 


The Euler solver uses a Ni-Lax-Wendroff node-based scheme on an unstructured grid. 
An explicit time-marching procedure subject to appropriate boundary conditions is used 
to drive the solution from an initial guess to a steady-state or to an unsteady solution. 

The numerical procedure has been introduced by Ni [47] for two dimensions and 
has been described later by Ni and Bogoian [48] for a three-dimensional application 
on a structured grid. The present chapter deals with the algorithm description for 
unstructured meshes. 

The spatial discretization uses hexahedral cells and the change in time of the state 
vector U is expressed as a function of the fluxes across the cell faces. These are evaluated 
as the average of the fluxes at the corner nodes. The residual (found from s ummin g 
the fluxes across the six faces of each cell) is used to determine the change in the state 
vector and is then distributed back to the nodes following the Lax-Wendroff algorithm. 

The state vector U at time-step n + 1 can be expanded in a Taylor series up to the 
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second- order terms as 


and the time derivative of the state vector U is related to the spatial derivatives of the 

fluxes F , G, H by Equation (2.1) 

dU _ _0F _ dG _ dH_ (3 

dt dx dy dz 

Hence, the change in the state vector between two iterations is 


.. rdF dG dH 1 

su = U n+1 -U n = - At — + — + — 

. ox dy dz . 

2 [dx V dt) 

The second-order changes are defined as 


d_ 

dy 




A F = At 


dF 

dt ’ 


AG = At 


dG 
dt ’ 


A H = At 


dH 
dt ' 


(3.4) 


Integrating over a pseudo- cell P formed by joining the centers of the cells surrounding 
node 1, as sketched in Figure 3.1a), gives 


J r W - l [-*< (f + f • + f )’ ■ - T {tx^ + Ty AG " + £**")] 

(3.5) 


Then by applying Gauss’ theorem we get 


SU 1 = -^ I (F,G,H)-ndS-^- f ( AF, AG, AH) • n dS. (3.6) 

Vj jp 2vx Jp 

n denotes the unit normal to the cell surfaces pointing outwards, V\ is the volume of 
P. Ati is the time-step associated with node 1 and is defined by 



(3.7) 


where the s um operates over the eight cells surrounding node 1. In Figure 3.1a) node 1 is 
surrounded by mesh cells A, B, G, D, E, F, G, H (mesh cells D and H are not represented 
for clarity purposes). Figure 3.1b) and c) represent the pseudo-cell split into eight cells 
Ap, Bp, Cp, Dp, Ep, Fp , Gp, Hp. The integration of the first-order terms is found by 
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b) 


B P 



c ) 

Figure 3.1: a) Mesh cells A to H surrounding node 1 and pseudo-cell P centered on 
node 1, b) and c) enlargement of pseudo-cell P split into eight cells Ap to Hp . Shaded 
surfaces indicate surfaces used for b) integration of first-order terms and c) integration 
of second-order terms. 
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integrating over all the surfaces of these eight volumes as shown by the dashed surfaces 
of Figure 3.1b) 


-TT- / (F,G,H)-ndS = -^p- / (F,G, H) • n dS. (3.8) 

V l J P Vi Ja p ,...,Hp 

Each of these integrals is estimated as one eighth of the surface integral over the corre- 
sponding mesh cell surrounding node 1 so that 


(F,G,H)-ndS = -^F [ (F,G,H) ■ fi dS. (3.9) 

V\ Jp 


The integral of the second- order terms is performed over the external surfaces of the 
eight cells Ap to Hp as represented by the dashed surfaces in Figure 3.1c), so that 

- ^ / (A P, A G, AH) ■ RdS = f (A F, AG, AH) ■ n dS. (3.10) 

2Vj Jp 2V\ J Api f Ap2iAp$ 

•ffp i > Bp2yHpy 


The integral over the surface Ap\ is evaluated as one fourth of the integral over 
the ‘mean surface’ ~A\ defined as the average surface between two opposite faces of mesh 
cell A. The ‘mean surface’ A 3 is represented on Figure 3 . 1 a). A similar procedure 
applies for surfaces Ap 2 to Hpz so that 

- t (A F, AG, AH) • HdS = [_ (AP, AG, AH) ■ fi dS. (3.11) 

2Vi Jp 8Ki JA u A 2 yAi 


Hence, 


s Ui = ~ f (F, G, H) ■ n dS - ^ L _ _ ( AP, AG, AH) ■ fi dS, (3.12) 
8 Vi J a H 8 Vi JA u A 7 yA> 




or formally 


SI 7i = SU 1A + SUib-*' + SU 1E , 


(3.13) 


where SU\a is the contribution from cell A to the change in the state vector at node 1 
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and so on for the contributions of cells B to H. SUi A cam then be written as 

su ™ = E_ ( af aS x + AG A S y + AJ?„S 2 )] , (3.14) 

\ A\ ,A2,A 3 / 

where Va is the volume of cell A. S x , S y and S z are the components of the surface 
vector in the Cartesian coordinates, and A Ua is the average first-order change in the 
state vector in cell A defined as 

A Ua = / {F,G,H)-ndS = -^ (FS X + GS y + HS Z ) (3.15) 

V A JA V A a 

and F, G , H are averages of F, G , H over the four nodes of each face. 

The second-order terms AF 4 , A Ga, A Ha are expressed as a function of the change AUa 
as 

AF ' = (m) A AU *’ ag *={w) a au *' ah * = {w) a av *- < 3 ' 16 > 

However, a more straight forward way to compute the second-orders terms is to use the 
changes in the conservative variables as follows. 


For a compressible flow 
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where ua, va, wa are obtained from averages over the nodes of cell A and 
(Au)yi = ((A(/m) - uAp)/p) A 

(At;),! = ({A(pv)-vAp)/p) A 

{Aw) a = {{A{pw)-wAp)/p) A 

(Ap)x = (7 - 1) ^A(pE) - uA(pu) - vA(pv) - wA(pw) + -^-(tx 2 + v * + w2 )^J 

For an incompressible flow 


AU a = 


’ A p* 


1 

n 

p w 

> 

e 

A u 

, A F a = 

2uAu + A p* 

Av 


uAv + vAu 

Aw 


uAw + wAu 


A 



r 



c 2 Av 


c 2 Aw 

vAu + uAv 

, A H a = 

wAu -|- uAto 
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AG a = 


where ua, va, tv a are obtained from averages over the nodes of cell A. 


In addition to Equation (3.14), the change in the state vector at node 1 receives 
contributions from the cells B to H written as 

6U lB = ^(^-AC/ b -_E_( A ^ 5 * + AG ® 5 v + A ^ 5 -))’ ( 3 * 17 ) 
1 V B B1.B2.B3 / 


1A ti I 
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SUic = tir-l#- AC ^- '£_(&FcS x + AG c S y + AHcS z ) , (3.18) 

1 \ c C1.C2.C3 / 

SU 1D = ^ (xr AC/ i>- E {&FdS x + AGpSy + AHdS z )\ , (3.19) 
1 V D di.dj.d, / 
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SU IE = Y, (&F E S x + AG E S y + AH E S Z ) ) , (3.20) 


E\ 1E2 }E$ 


SU 1F = Y ( af fS x + AG F Sy + AH f S 2 ) ) , (3.21) 


F 1 ,^2,^3 


SU lG = r 


1A«j / V G 


8 Vi \ At G 


A U G - Y (& f gS x + AG G S y + AH g S z ) , (3.22) 


Gi,G 2 ,G, 


SVlB = T 


1 At, / V H 


8 V] l Atjj 


A/7* - ]T (AFjjS x + AGffSy + A H n S z ) ) . (3.23) 

HuH 2 ,Hy 


The calculation of the cell volumes, face areas and volumes associated with cell nodes 
is described in Appendix B. 
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3.2 Numerical smoothing 


By expanding the dispersion relation of the Lax-Wendroff scheme into Taylor series, it 
can be seen that the Lax-Wendxoff scheme carries both an inherent dissipation term 
and a third-order dispersion error [73, 76]. The former term accounts for the scheme 
stability when operating on smooth flow fields. Nevertheless, the dispersion error is 
responsible for the introduction of oscillatory modes in the solution. In order to damp 
these background oscillations, an artificial dissipation term has to be added to the Euler 
equations. 

The compressible flow version of the Euler solver uses a standard Laplacian second- 
difference smoothing. The Laplacian is obtained at node 1 by su mm i n g the difference 
between the state vector at node 1, I7i, and the cell-averaged state vector of the 8 
surrounding cells Ua to Ujj (the nomenclature is described in Section 3.1). The change 
in the state vector at node 1 is then the sum of Equation (3.13) and the contribution 
from the second- difference smoothing 

+ + + E (3-24) 

where is an artificial viscosity coefficient, V\ the volume associated with node 1 and 
Vkc the volume of the cell kc. Ati and Ai*. c are the time-steps associated with the node 
1 and the cell fcc, respectively. 

The incompressible flow version of the Lax-Wendroff algorithm uses a fourth- difference 
smoothing. The fourth- difference operator is constructed as a second- difference of a 
second-difference. Instead of using two standard Laplacian operators, the inner second- 
difference involves a “pseudo-Laplacian” operator defined by Holmes and Connell [31] 
and obtained at node n as 

tmfli 

D 2n = E MUi - U n ), (3.25) 

»=1 

where the sum operates on any number i max of nodal points surrounding node n. As 
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Figure 3.2: Nodal points 1 to 6 chosen for determination of fourth- difference smoothing 
at node n. 


shown in Figure 3.2 the 6 surrounding nodes along the neighboring cell edges have been 
selected in the domain (less nodes are involved at the boundaries because of the missing 
neighbors). The Wi are the “weight” factors allowing the smoothing operator to be 
transparent when applied to any linear function in x, y or z and for any choice of i max . 
Thus, the weight factors must statisfy 


1-max 

Yi w '( x ' ~ *«) = 

*=i 

(3.26) 

*mdi 

Z2 W i(Vi - yn) = 0, 
1=1 

(3.27) 

*mat 

2Z W i( Z * ~ Z n) = 0. 

i-l 

(3.28) 


To ensure stability, the weight factors are chosen close to unity by minimizing the cost 
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function 


*m.ax 

C n = & w h ( 3 - 29 ) 

i=l 

where 

A W{ = 1 — Wi. (3.30) 

The fourth- difference smoothing is written as the standard second- difference operating 
on this “pseudo-Laplacian” . The change in the state vector at node 1 is then formed by 
the sum of Equation (3.13) and the contribution from the fourth- difference smoothing 

6U\ = SUia + ••• + SU\h - V *~7T~ 'TT~(^ , 2 kc ~ D 2i )- (3.31) 

Vl kc=A,...H alkc 

i >4 is the fourth- difference smoothing coefficient and D 2kc is the cell-average of D 2 at 
the eight nodes of cell kc. 

The use of this pseudo-Laplacian ensures that the smoothing term does not distort 
any linear function in the domain (even the linear functions perpendicular to the bound- 
aries). Since the choice of i max is arbitrary, there is no need for a special treatment of 
the smoothing at boundaries of the domain (inlet, exit, wall,...). Also, the use of weights 
allows linear functions to be exactly represented on irregular grids. 

In some instances, as in the case of a convex wall, some of the weights become 
negative. Holmes [31] elected to clip the weights in the range (0,2) because of stability 
problems. However, it has been found in this work that clipping the weights resulted in 
a degradation of the solution. The negative weights have not been clipped in this work 
without any sig nifi cant stability problem. However, a slight decrease in the convergence 
rate was observed. 

The actual algebra involved in the computation of the weight factors is presented 
below. The problem consists in minimi zing the cost function C n under the constraint 
that the pseudo-Laplacian operator must give zero value when applied to linear functions 
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in x, y and z. Using the Lagrange multipliers, the function to minimize is 

imax 

^ ~ W i{^x{ x i ~ x n) + A y (t/* ~ Vn) + A z (z^ ~ z n))) 


i = 1 


(3.32) 


where A x , A y , A,, are three unknowns. The ‘Euler’ equation of the problem is 

dF 


dAw{ 


-0, 


or 


Awi — 0.5(\ x (xi — £ n ) + A y (y{ — y n ) -f A Z (z{ — z n )) 


(3.33) 


(3.34) 


By combining this equation with Equation (3.26),Eq (3.27) and Equation (3.28), A x , A y , A, 
are found as 
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(3.35) 
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3.3 Mesh singularity treatment 


The meshes used in this work for the computation of flow in pipes present singularities 
where a given node belongs to only six cells instead of eight. Figure 3.3 represents such 
a singularity (in a 2-D case for clarity). Node 1 is surrounded by cells A to F (faces of 
cells A to C are shown). The pseudo-cell P has volume Vi defined as an average over 6 
cells 


Vi = J £ V all . 

6 6 cells 

Ati is the time step associated with node 1 defined by 



(3.36) 


(3.37) 


The pseudo-cell P can be split into 6 smaller cells (faces of Ap to Cp are shown). 



Figure 3.3: Mesh singularity at node 1 with pseudo-cell P. 
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The integration of the first-order terms is performed over all the surfaces of these 6 
cells. As in Section 3.1, each of these integrals is in turn estimated as one eighth of the 
surface integral over the corresponding mesh cell surrounding node 1. Therefore 

f (F,G,H)-ndS = -^- [ ( F,G,H)-ndS . (3.38) 

The integral of the second-order terms is performed on the external surfaces of the 6 
cells Ap to Fp (surfaces Api, Ap 2 , Ap 2 to Fp\, Fp 2 , Fp 2 ). The integral over Ap\ is 
evaluated as one fourth of the integral over A\ and so on for faces Ap 2 to Fp$. Hence, 


- ^ 7 - f (AF,AG,AH)-ndS = [_ (AF,AG,AH)-ndS. ( 3 . 39 ) 

J P 0V1 JA x ,A 2 i A 3 


Fit F 2 , Fy 


According to Section 3.2 the additive smoothing term can now be written as a sum over 
6 cells instead of 8 and the final change in the state vector at the singular node 1 is 


Ati 


SU\ — + 6U\b + 6Uic + &U\d + SUie + + ^2-fr- 'Y' ( Y kc {Uk c - 

V l 


kc=A,... y F 


Ui))> 

(3.40) 


when using a second-difference smoothing and 


SUi — SUia + SUib + SUic + + SUie + SUif - 1/^—^- V] {D 2k - d 2i ), 

V 1 kc=A,...F Atk ‘ 


(3.41) 

when using a fourth- difference smoothing. 6U\a to &U\F are given by Equations (3.14) 
and (3.17) to (3.21). 


In the unstructured code implementation the update of the state vector for the sin- 
gular node is computed in a loop over the cells. Hence, the fact that only six cells are 
contributing to a singular node is directly dealt with by the connectivity table. Never- 
theless, V] and At-[ associated with the singular node 1 need to be defined according to 
Equations (3.36) and (3.37). 
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3.4 Farfield boundary conditions 


The farfield boundary conditions allow the use of truncated domains without affecting 
the n um erical solution. These conditions follow from the 1-D linear characteristics 
theory whose theoretical background and implementation are described by Giles in 
[23] and [22]. The derivation of the Euler equations written in primitive form and 
computational coordinates is described in Appendix C. These equations are readily 
transformed back to physical coordinates as 


dt + dx 


+ B 


8y dz 


= 0 . 


(3.42) 



u p 0 0 

0 u 0 0 

A = 0 0 u 0 

0 0 0 u 

i 0 7p 0 0 






cl 0 
2u 0 
V u 
w 0 
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/ o 0 cl 0 \ 
0 v u 0 


B = 


1 0 2v 0 

y 0 0 w v J 

for an incompressible flow (the artificial 
the Euler equations). 


^000 
0 w 0 

, c = 

0 0 w 

^10 0 

compressibility method 



u 


v 

2w ) 

has been used to modify 


Considering perturbations from a uniform and steady flow and retaining only the 
first-order terms gives a linear equation for J7 p , the vector of the perturbations from a 
uniform flow 


°Z> + a M’ +b Ml + c °&- 


= 0 , 


(3.43) 


dt dx dy ' ~ dz 

where .4, B and C are based on the uniform steady flow. Furthermore, if it is assumed 
that the perturbations travel normal to the considered boundary (say in the x direction), 
then the derivatives in the y and z directions can be ignored and the above system of 
equations reduces to the 1-D linear system 




dt 


dx 


(3.44) 


For convenience, the* symbols will be omitted in the remainder of this section and U p is 
renamed t/. The convention is that the matrices 4., B and C are based on the uniform 
flow field and that U is now the vector of the perturbations from the unif orm flow field. 
The matrix A can be diagonalized by the transformation 


A = T~ l A T, (3.45) 

where the columns of T are the right eigenvectors and the lines of T -1 are the left 
eigenvectors associated with the five eigenvalues of A. Multiplying Equation (3.44) by 
T -1 gives 

T" 1 ^ + T- 1 ATT- 1 ^ = 0, (3.46) 
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or 


d$ n 

li + A dx~ ’ 


(3.47) 


where $ is the vector of the linearized characteristic variables defined as 


$ = T~ l U. 


(3.48) 


The eigenvalues of A (i.e. the diagonal elements of A) represent the speeds of different 
propagating waves. Each positive eigenvalue corresponds to a characteristic entering the 
domain and each negative eigenvalue corresponds to an outgoing characteristic. At any 
boundary the number of conditions to be imposed must correspond to the number of 
characteristics entering the domain. 

In the following the superscript * stands for specified values. Superscript n refers to 
values at the previous time-step and n+1 to predicted Lax-Wendxoff values at the present 
time-step. Any characteristic leaving the domain will be superscripted n+1 since it is a 
function of the Lax-Wendroff values in the domain. 

3.4.1 Farfield boundary conditions for compressible flow 

In the case of a compressible flow, the matrices T and T -1 are 


1 

0 

0 

1 

2? 

1 \ 
2? 


( -c 2 

0 

0 

0 


0 

0 

0 

1 

2 pc 

1 

2 pc 


0 

0 

pc 

0 

0 

0 

1 

pc 

0 

0 

0 

, T- a = 

0 

0 

0 

pc 

0 

0 

0 

1 

pc 

0 

0 


0 

pc 

0 

0 

1 

V 0 

0 

0 

1 

2 

1 

2 ) 


V 0 

-pc 

0 

0 



and the five eigenvalues of matrix A are 

Ai = u, A 2 = u, A 3 = u, A 4 = ti + c, A s = u - c. (3.49) 
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For a subsonic compressible flow the first four eigenvalues are positive and the fifth 
is negative. Thus, at an inlet boundary the four first characteristics are propagating 
into the domain whereas the fifth propagates out of the domain. At an exit boundary 
the first four characteristics are exiting the domain and the fifth propagates into the 
domain. 


Inlet boundary - compressible flow 


At the inlet boundary, the total enthalpy /z 0 , the entropy related function s and two 
angles ai, c*2 of the inlet velocity vector are specified. This is achieved by defining the 
four residuals Ri to R4 as differences between specified and Lax-Wendroff values as 


R\ — — ^0’ 

R2 = s n -s 5 , 

iZ 3 = (tana x ) n — (tan ax) 3 , 

R4 = (tana2) n — (tana^)*, 

The entropy-related function s and the angles a x ,a2 are defined as 


s = log(7p) - 7 log p, tan a x 


v w 

— , tana 2 = — . 
u u 


The necessary changes in the characteristics in order to drive the residuals to zero are 
found by performing one step of a Newton-Raphson procedure as 


f Ri > 


f Si 1 ) 

R 2 

d(Rj, R 2 , R$, R 4 ) 

d($a,$ 2 ,^3,$4) 

Si 2 

£3 

Si 3 

R 4 j 


^ ^§4 y 


The Jacobian matrix is computed as 


( 3 . 50 ) 


J = 


d(Riy R 2 y # 4 ) _ d[Ru R 2 , R3, R*) d(p,UyVyW,p) 

d($!, $2, $3, $4) d(pyUyVyWyp) ^3,^4) 
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V w 1 c + u 

pc pc pc 2 

0 0 0 

n -v 
P cu 2 pcu 2 

n 1 -tt> 

P cu 2^? 



The required changes in the first four characteristics are then found from Equation (3.50) 



whereas the change in the fifth characteristic is found from 


(3.51) 


(£$ 5 ) n+1 = Tf 1 SU = —pc(Su) n+1 + {6p) n+ \ (3.52) 


where Tf 1 is the fifth line of T _1 and (<5p) n+1 and (rfu) n+1 are the Lax-Wendroff changes 
in pressure and inlet normal velocity at the inlet nodes. The changes in the primitive 
variables are then found by 
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Exit boundary - compressible flow 


At the exit boundary the pressure is set corresponding to the one necessary condition 
defining the residual 


R = p n -p s . (3.53) 

The necessary changes in the characteristics in order to drive this residual to zero is 
found by 

dR 

r + w ; s * 5 = 0 ’ (3 - 54) 

and 

OR _ dR d(p, u, u, w,p) 1 
d$ s d{p,u,v,w,p) d $ 5 ” 2 ' 

The required change in the fifth characteristic is then found as 


= -2Jii = -2 (p n -p 5 ), 


whereas the change in the first four characteristics are calculated using the Lax-WendrofF 
changes in the primitive state vector as 


f 6$ i 

S$ 2 

6* 3 


n+ 1 


V ) 


( -c 2 

0 

0 

0 


0 0 0 

0 pc 0 

0 0 pc 

pc 0 0 


1 

0 

0 

1 


\ 



6u 


6v 


Sw 

/ 

Sp / 


n+1 


The changes in the primitive variables are then found by 


6p > 


1 (£$ 1 ) n+1 ^ 

6u 


(** 2 ) n+1 

6 v 

= T 

(*$ 3 ) n+1 

Sw 


(6$ 4 ) n+1 

6p y 


^ ^$5 j 


(3.55) 


51 



The changes in conservative variables are found for an inlet or exit boundary by 
using the changes in the primitive variables as 

S(pu) = u6p + pSu, 

S(pv) = vSpi-pSv, 

6{pw) = wSp + pSw, 

S(pe 0 ) = — - — 6p + puSu + pvSv + pwSw + -(u 2 + v 2 + w 2 )Sp. 

7 — 1 * 


3.4.2 Far field boundary conditions for incompressible flow 


In the case of an incompressible flow, the matrices T and T 1 are 

/ 


/ 


T = 


0 0 


c 2 J 3 


-cl(3 


\ 


0 0 u(u + /3) + cl u(u-/3) + cl 

0 1 v(u + /3) v(u- P) 

^10 w(u + P) w(u-P) ) 


T ~ 1 = 


w 

p 2 

V 

'W 


wu 

p 2 

vu 
2 


P 

t j-P 1 
2 p' 

u+J. 3 _L 

V 2/3 


and the four eigenvalues of matrix A are 

Aj = u, A2 = tt, A3 = u + Pi A4 = u — P, 


0 1 ^ 
1 0 
0 0 
0 0 


(3.56) 


where P = \/u 2 + c\. 

The first three eigenvalues axe positive and the fourth is negative. Thus, at an inlet 
boundary the first three characteristics sire propagating into the domain whereas the 
fourth propagates out of the domain. At an exit boundary the first three characteristics 
are exiting the domain and the fourth propagates into the domain. 
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Inlet boundary - incompressible flow (internal flow case) 


For the internal incompressible flow cases, namely the flow in pipes, the inlet velocity 
profile is known through measurements or can be approximated by computing a fully 
developed velocity profile. In these cases, the three components of the inlet velocity are 
specified so that 


f?i = u n — u s 
R 2 = v n - v s 
Rz = w n — w 3 


The necessary changes in the characteristics in order to drive the residuals to zero are 
found by 


: 


d(-Rl, -R2; -R3) 




^2 

j 


= 0. 


The Jacobian matrix is computed as 


(3.57) 


j _ d(-Ri, #2,-^3) _ d(p*,u, u,u?) 

#($i, $2, $3) d(p*, u , v, w) d($i, $ 2 , $3) 


/ 


\ 


0 1 
0 0 
0 0 


0 0 ^ 
1 0 
0 1, 


{ 0 0 
0 0 
0 1 
K 1 0 


cl 3 > 

u(u + 0 ) + cl 
u(u+ /3) 
w(u + (3) j 


^00 u(tx + /3 ) + cl ^ 
0 1 v(u + (3) 

^10 w(u + /3) j 


The required changes in the three first characteristics are then found from Equa- 
tion (3.57) as 


( 

6$ 2 
y ^$3 


\ 


/ 


-J 


-1 


( Ri \ 

R 2 

\ R * } 


(3.58) 
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whereas the change in the fourth characteristic is found from 


(**4) n+1 = T. f 1 SU = -^±P(6p*)”+' + - ^(Su) n+ \ (3.59) 

77 1 is the fourth line of T~ l and (6p*) nJrl and (£u) n+1 are the Lax-Wendroff changes 
in pressure and inlet normal velocity at the inlet nodes. The changes in the primitive 
variables are then found from 


Sp* ^ 


( } 

6u 


6$ 2 


= T 


6v 


^3 

Sw j 


V (^4) n+1 ) 


Inlet boundary - incompressible flow (external flow case) 

For the external incompressible flow cases, namely the flow around a wing, better results 
were obtained by not specifying the three components of the inlet velocity but the 
stagnation pressure and two angles of the velocity vector instead. In this case the 
residuals are 


R\ = Po “ P01 

R 2 = (tanai) n — (tanai)% 

R z = (tana 2 ) n - (tana 2 )% 

where the stagnation pressure and the two angles of the inlet velocity vector are defined 
as 

* 1, 2 2 2x v w 

Po = P + ~(u z + v l + nr), tanai = — , tana 2 = — . 

2 u u 
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The necessary changes in the characteristics in order to drive the residuals to zero are 
found by 


R 2 

U» / 


djRiy R2, Rz ) 

0 (*1,*2,* 3 ) 


/ 


\ 


8$ 2 
8 $ 3 ) 


= 0. 


The Jacobian matrix is computed as 


J _ d(Ri, i? 2 , ) _ d(Rj, i?2; -R 3 ) #(?*, v, w) 

5($1, $2, $3) d(p*, U, V, w) d($i, $2, $3) 


(3.60) 


/ 


\ 


1 

0 

0 


u 


tan Qt 

u 

tan 0:2 

u 


v 

i 

u 

0 


A J 


f 0 0 
0 0 
0 1 
\i 0 


c\P ' 

u(u + P) + cl 
v(u + P) 
w{u + P) 1 


w t? (u + P){c 2 a + u 2 (l + tan 2 Qi + tan 2 a 2 )) \ 


0 i 
v i 0 


tanajc^ 

u 

tana^ 

u 


/ 


The required changes in the first three characteristics are then found from Equa- 
tion (3.60) as 


( £$1 ^ 


8* 2 
8$ 3 


= -J- 1 


1 

R 2 

\ R 3 / 


(3.61) 


whereas the change in the fourth characteristic is found from 


(** 4) n+1 = T ^ 8U = -^|(*PT +1 + 2^(^)" +1 . (3.62) 

T" 1 is the fourth line of T~ l and (6p*) n + l and (6u) n+1 are the Lax-Wendroff changes 
in pressure and inlet normal velocity at the inlet nodes. The changes in the primitive 
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variables are then found from 


6p* ^ 



6u 

= T 

6*2 

6 v 


6$ 3 

6w j 


, (** 4 )" +1 j 


Exit boundary - incompressible flow 


In the pipe flow cases, the vortex created by the secondary flow impinges on the exit 
boundary* Also in the incompressible flow around a wing, the trailing vortex sheet rolls 
up into a vortex which crosses the exit boundary. Imposing a constant pressure at the 
exit boundary would not be consistent with the presence of the vortex. Instead, good 
results were obtained by setting the first derivative of pressure to zero in the direction 
normal to the exit boundary. Using the node i — 1 adjacent to the exit boundary as 
depicted in Figure 3.4, this gives 

_* 

Pi = Pi-i 

Approximating the pressure in the direction perpendicular to the exit boundary by a 
linear function (by setting the second derivative of the pressure to zero) gave similar 
results. 

Using the pressures computed at each node of the exit boundary as specified pres- 
sures, the following residual is defined 

R=p* n - p“. (3.63) 

The necessary changes in the characteristics in order to drive this residual to zero is 
found from 

7 ? 

R + = °> ( 3 - 64 ) 

09 4 
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vortex axis p ~ cst 



z — 1 

z exit boundary 


Figure 3.4: Exit boundary with impinging vortex. 


and 


dR 


dR d(p*,u,v,w) 


= 1 ■{-&). 


<9(p*,u, v,w) 

The required change in the fourth characteristic is then found as 

^§ 4 = = - 4 ( p *" - p **), 

cl(3 cl/3 yF P h 

whereas the change in the first three characteristics are 


(^$ 1 ) n+1 = Tf 1 SU 

_ w 


(fa )** 1 

(S$ 2 ) n+1 = Tf 1 SU 

II 

4* 

(<pT +1 - p 


(** 3 ) n+1 = T 3- 1 SU 

_ -«+ / 3 

(W + 5P 

(*u) n+1 . 

2/3 2 c 2 




n+ 1 


Xf 1 , T^" 1 , X^ 1 are the first, second and third lines of X -1 , respectively and (<5p*) n+1 , 
{Su) n +\ (£t;) n+1 , (^tt7) n+1 are the predicted Lax-Wendroff changes in pressure over den- 
sity and velocity components at the exit nodes. The changes in the primitive variables 
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are then found from 


8p* ^ 


( (^$ 1 ) n+1 ^ 

6u 

= T 

(£§ 2 ) n+1 

8v 


(** 3 ) n+1 

8w j 


y ^$4 j 
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3.5 Wall boundary condition 


A no flux boundary condition is imposed at the wall boundaries, i.e. 

v *n = 0. (3.65) 

This condition implies that the fluxes F^G , H through the wall faces are only a function 
of the pressure. 


The wall boundary condition is imposed once the Lax-Wendroff fluxes have been 
computed and the changes in the state vector have been determined for each node. 
Since the boundary condition deals with a flux condition on a face, the cell residual 
computed before in the main Lax-Wendroff routine is erroneous. The wall boundary 
condition corrects the residual for these particular cells. The corrected cell residual is 
then distributed to the eight nodes as before. The cell residual for cell A is (taking into 
account only the first-order terms) 

A U A = -y± J2 (■ FS x + GS y + HS z ), (3.66) 

^ 6 faces 

where F, G, H stands for the average values over the four nodes of F, G, H. The residual 
at the Weill node 1 of cell A is defined as 


- 5TT (^ Va ) ■ (3 ' 67) 

For a wall node, the node residual receives contribution from only four cells instead of 
eight but also takes into account a smaller node volume. By adding to the cell residual 
A V a the term 



0 ) 
pS x 


FS X + GS y + HS Z - 

pSy 

ps z 



\ 0 ) 
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computed on the wall face for the compressible flow cases and the term 

A t_A 
V A 

computed on the wall face for the incompressible flow cases, the contribution of a normal 
mass flux on the wall face is eliminated. The new cell residual is then distributed to the 
eight nodes as in Section 3.1. 



0 1 


FS X + GSy + H S z - 

P*S X 

¥Sy 



V p 5 ^ ) 
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3.6 Symmetry boundary condition 


In the incompressible flow cases analyzed in this work, the symmetry of the problem 
allows the solution to be performed on only half of the domain. A symmetry-plane 
boundary condition is enforced on the nodes lying in the symmetry surface. 


The symmetry surface is defined by one computational plane. Due to the present 
choice of Cartesian coordinates, this surface corresponds also to the (t/ = 0.) physical 
plane. Then for any existing node at a position ( x , y, z ), a pseudo-node located at a 
symmetrical position ( x , —y, z ) is characterized by the same primitive variables but an 
opposite velocity component in the direction perpendicular to the symmetry surface 
(here the v component). For a given face with vector surface 5 = (5 x ,S y , 5 Z ), the 
surface vector of the symmetrical face is 5, = (S x , -S y , 5 Z ). From Equation (3.15) the 
average first-order change in the state vector A U a , for the pseudo-cell A, symmetrical 
to cell A is then 


' A p* ^ 


Ap* \ 

A u 


Au 

Av 


-Av 

^ Aw j 

A. 

Aw / 


From Section 3.1, it can be seen that the same kind of relations between cell A and cell 
A a holds for the second-order terms. The time-steps and the volumes of cell A and A s 
are identical, so that the contribution from cell A s to the change in state vector at node 
1 lying on the symmetry surface is 


(s P - \ 


' 6p * ^ 

6u 


Su 

Sv 


—Sv 

y 6w j 

1A § 

6w / 


Therefore, the contribution from symmetric pseudo-cells is taken into account by 
doubling the contribution from the existing cells for the pressure and the components 


61 



of velocity parallel to the symmetry surface. The change to the component of velocity 
perpendicular to the symmetry surface is set to zero. 

Because the initialization of the flow may not satisfy the symmetry condition, a ‘no 
through velocity’ condition is specifically applied each time the state vector is updated, 
that is vi = 0. 

The smoothing stencil at the symmetry surface is completed by taking into ac- 
count the pseudo-node at location (z, — t/,z) characterized by the primitive variables 
p*, u, — v, w. 
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3.7 Numerical implementation on unstructured grids 


The actual implementation of the Lax-Wendroff algorithm is done in two steps. First 
the fluxes F,G,H at the nodes are computed in a loop over the nodes. Secondly, in 
a loop over the cells, the cell average state vector is computed. In the same loop the 
changes due to the first-order and second-order terms are determined by summation over 
cell faces, and the first-order and second-order terms are then added to the smoothing 
term to give the contribution of each cell. A second- difference operator is computed in 
the same loop and is applied on the state vector U in the case of a second-difference 
smoothing operator. In the case of a fourth difference smoothing, this operator is applied 
on a pseudo-Laplacian computed beforehand in a separate routine. The change in the 
state vector for the 8 nodes is then obtained by summing the contributions from the 
eight surrounding cells. 

The subsequent boundary conditions imposition modifies the changes in the state 
vector for the nodes lying on boundary surfaces. Finally, in a loop over the nodes, the 
state vector changes are added to the current state vector value for each node. 

The Lax-Wendroff algorithm is implemented on an unstructured grid where an array 
over the cells points to the 8 nodes of each cell. When looping over the cells, this array 
allows one to readily add the contributions from each cell to a node. 6 arrays over 
the cells are initialized as the 6 surface vectors for each cell. The numbering of the 
surface vectors is described in Figure B.l of Appendix B. The inlet, exit, wall and 
symmetry nodes have associated arrays to allow the differentiation of these particular 
nodes from the rest of the field when handling the boundary conditions. The fourth- 
difference smoothing term calculation requires ‘neighbor 5 information. According to 
the stencil defined in Figure 3.2, an additional array over the nodes contains the 6 
neighbouring nodes along cell edges. Another array contains the 6 weight factors used 
in the computation of the pseudo-Laplacian. 

During the calculation each node gets associated with the fluxes F, G, #, a volume, a 
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time-step, a change in the state vector and the state vector itself. Each cell is associated 
with a volume and a time-step. 

In order to take advantage of vectorization capabilities, each cell is ‘colored’. The 
loops over cells are replaced by loops over colors. The ma x i m u m number of colors is 
determined by the requirement that any two cells of one color must have no common 
nodes, so that dependencies can be avoided. Boundary cells are associated with separate 
colors in order to differentiate them from the field cells. 
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3.8 Time-step restriction 


According to the stability analysis given in Appendix C, the allowable time-step is given 
by the CFL condition 


At < 


yjr j + r| + r| + cV a 2 + 6 2 + d 2 


for a compressible flow and 
At < 


\J T l + r 2 + r 3 + Vm ' + r 2 + r 3 + c a(° 2 + & 2 + d 2 ) 
for an incompressible flow, where 


(3.68) 


(3.69) 


J = Xf(y n z c - y c Zr,) - x v (i f t Z( - y ( Z() + x ( (y (Zv - y v z ( ), (3.70) 

n = u(y„z c - y <Zv ) - v{x vZ<: - x (Zv ) + w(x v y ( - x ( y v ), 

ri - -t i(y ( z c - y (Z( ) + v(x (Z< - x (Z( ) - w{x^y <: - x^), 

r 3 = u{y tZr) - y„z ( ) - v(x iZv - x vZc ) + 10(3;^ - x v y ( ), 

« 2 = {y^c-y^nf + iy^c-ycHf + iyc^-yr,^) 2 , (3.7i) 

b 2 = (x vZ( - x (Zr ,) 2 + (x (Z< - X (Z( ) 2 + (x (Zr] - x v z ( ) 2 , (3.72) 

^ = i x r ] yc-x <i y i ) 2 -\-{x i y (i -x < y t ) 2 ^(x i y r ,-x v y t ) 2 . (3.73) 


The CFL condition is satisfied with a margin set by the CFL number multiplying 
the ma x i m um time-step. This number is set to a value of 0.9. When computing steady- 
state flow problems with the time-marching procedure, local time-stepping is used in 
order to accelerate the convergence, i.e. the CFL condition is satisfied independently 
for each cell in the flow field. When computing unsteady flow cases, the time-step used 
in each cell is set equal to the mini mum time-step in the field. 
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3.9 Accuracy study 


An accuracy study is performed to confirm the second-order accuracy of the Eule- 
rian Lax-Wendroff scheme. Since the incompressible flow through bent pipes is one of 
the topics of this thesis, the constant stagnation pressure flow though a 90° bend is 
chosen for this accuracy study. The Euler equations are transformed using the pseudo- 
compressibility concept described in Chapter 2 and since the flow is symmetrical, the 
calculations are performed on only one half of the pipe. 

The test geometry is the 90° bend of circular cross-section, taken from the Enayet 
et 2 d. data set [21]. It will also be used in Chapters 10 and 11 of the results. The pipe 
diameter is 0.048 m and the ratio of radius of curvature to pipe diameter is 2.8. The 
geometry extends two diameters upstream and downstream of the bend. The accuracy 
study is performed by using three grid densities with a mixed O-H type grid on each 
pipe cross-section and a H type grid along the bend. The three grids are composed of 
53 x 22 nodes, 189 X 43 nodes and 713 X 85 nodes, respectively. Each grid is composed 
of four times as many cells as the previous grid in a cross-section and two times as many 
cells in the streeimwise direction. Side 2 ind front views of the grids are represented in 
Figure 3.5. 

The inl et conditions 2 ire defined to have constant stagnation pressure. A constant 
inlet velocity is set through the inlet boundary condition. The inlet cross-section is 
placed two diameters upstream of the bend so that the upstreeim influence of the bend 
is negligible at this location and the constant inlet pressure condition is obtained. 


The exact solution for the flow field has constant stagnation pressure everywhere. 
Therefore, a globed error in the Eulerian solver is quantified in terms of the £2 norm of 
the errors in stagnation pressure e Po 



where N is the total number of grid nodes, and A is the local error in stagnation 
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Figure 3.5: Front and side views of the grids used in the accuracy study (53 X 22 nodes, 
189 x 43 nodes and 713 X 85 nodes). 
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Figure 3.6: L 2 norm of errors in stagnation pressure as a function of the grid spacing 
for the three grid densities. 


pressure coefficient defined at a node by 


AC P0 



(3.75) 


The denominator represents the inlet dynamic head and po in is the inlet stagnation 
pressure. Referencing to the inlet dynamic head instead of the inlet stagnation pres- 
sure eliminates the dependency on the background pressure level which is arbitrary for 
incompressible flows. 


The results are represented in Figure 3.6 where the L 2 norm of the errors in stagna- 
tion pressure coefficient is plotted as a function of the grid spacing h in a logarithmic 
scale for the three grids. The grid spacing is not uniform for this test case, but since 
the value of h is just a reference value, it is determined as an average over the cells. 
A line is drawn through the three points using a least-square fit. The slope of the line 
deteiniining the order of accuracy of the Lax-Wendroff scheme is -2.0. The theoretical 
value which ensures second-order accuracy is -2. 
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Chapter 4 


Lagrangian Governing Equations 


4.1 Lagrange equations for compressible flow 


The Lagrangian equations for the convection of a particle in an inviscid compressible 
and isentropic flow subjected to a conservative force field are 


DU i 
Dt 


= T U 


( 4 . 1 ) 


where 


Si s l + ('- v > < 4 - 2 > 

defines the material derivative. When applied to a fluid quantity, this operator expresses 
the ‘convective change’ of the quantity due to the displacement of a given particle in 
time as the sum of the unsteady and spatial changes as the particle moves from one 
location to another. Ui is the state vector of the Lagrangian variables and T\ is a source- 
term vector governing the change during time of the Lagrangian state vector for a given 
particle. The Lagrangian state vector Ui and the source term Ti are defined by 



r 


V 

Ui = 

2 

, T t = 

(<I5 • V)tT- ul(V • v ) - V(i) x Vp 


S 


0 
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where r — (:c,y, z) is the position vector. The vorticity u? and the entropy-related 
function S axe defined by 

(2 — V X v , S — — . (4-4) 

P 7 

By using the continuity equation, the change of vorticity for a given particle during time 
can also be written as 

“5^ = (2 . V)<r- iv(i) X Vp. (4.5) 

The Lagrangian view of the momentum and energy equations given by Equations (4.1) 
and (4.3) expresses the convective change of vorticity and entropy for a particle. 


4.1.1 Source-term contribution in compressible flow 

The convective change for in a compressible flow is the sum of the contributions 
from a tilting/stretching term and a non-barotropic term. The first term (^ ■ V)v results 
from the tilting and the stret chin g of the vortex lines during time. The non-barotropic 
term -iv(^) x Vp represents the contribution from the moment of the pressure forces 
about the center of mass. The entropy-related function 5 is said to be passive, since it 
does not undergo any convective change. 


4.2 Lagrange equations for incompressible flow 


For an incompressible flow subjected to a conservative flow field, C/j and Ti are defined 

by 



■ 


- 1 

Ui = 

r 

lS 

m - 

, T,= 

p > 
1 


(4.6) 
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The use of the total pressure po as an additional convected quantity has not been elected 
since the vorticity and the Eulerian velocity field already determine the total pressure 
gradient according to Crocco’s equation for incompressible flow 


dv 

It 


- V X u> = 



(4.7) 


4.2.1 Source-term contribution in incompressible flow 


In the case of an incompressible flow, the source-term reduces to the tilting/stretching 
contribution (u; • V)tT. By decomposing (u> • V)u, the stretc hin g terms are identified as 
the terms 


du dv 



(4.8) 


The fluid elements in the z, y, z directions with associated vorticity components u! x ,u? y ,u? z , 
are stretched due to the strain fields Qjj, respectively. The stretching of a fluid 

element of length dy with associated vorticity w y under the strain ^ is illustrated in 
Figure 4.1. The resulting change in vorticity is ( u? y dv ) for a fluid element of length dy. 
Each component of w is intensified or reduced depending on the the sign of the strain 
field. A positive strain corresponds to a vortex line intensification and inversely. The 
remainder of the terms 


du 


du 


dv 


dv 


u?. 


y dy ’ * z dz " X dx Wz dz 








dw 




dw 

dy 


(4.9) 


are due to tilting. An example of vortex line tilting is illustrated in Figure 4.1 where a 
fluid element of length dx with associated vorticity u> x is tilted an amount dw about the 
y axis under the influence of the strain field dw/dx . The resulting change in vorticity 
for a fluid element of length dx is w x dw. 


As can be seen in Figure 4.1, the final directions of the fluid element and vorticity 
vector are identical. Therefore, in an inviscid incompressible fluid subjected to a con- 
servative force field, the vorticity Q is said to be frozen to the fluid element. The fact 
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z 



Figure 4.1: Stretching of vortex line in y direction and tilting of vortex line in x direction. 


that vortex lines move as material lines, in an incompressible and inviscid flow, can also 
be deduced from the comparison of the Helmholtz equation 


D<2 _ 

— — = u> • Vv, 

Dt 

to the equation for the evolution of an infinitesimal line element dl [8] 

Ddl 


(4.10) 


Dt 


— dl • Vv. 


(4.11) 


In a compressible flow, the quantity — is attached to the fluid element. 
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Chapter 5 


Eulerian/Lagrangian Integration 


The Lagrangian solution is based on particle markers convecting through the flow. 
Each marker position, vorticity (and additionally entropy for compressible flow) at the 
marker’s position form the Lagrangian state vector [/). The convective change in the 
Lagrangian state vector, as the marker moves through the Eulerian grid, is obtained by 
the evaluation of the source term Ti. 

The Eulerian solution provides the information required for the computation of the 
markers’ trajectories and for the integration of the Lagrangian source-terms. In turn 
the Lagrangian solution is used to correct the Eulerian solution and reduce its numerical 
diffusion error. This process allows the Lagrangian solution to accurately capture the 
convection of vorticity (and entropy) while the Eulerian solution conserves mass, mo- 
mentum (and energy). The standard time-integration of the Eulerian solver supports 
the iterative interaction of the Eulerian solution and the Lagrangian particle tracking 
solution. 

First the computation of the markers’ convection and the integration of the source- 
terms are described. The correction of the Eulerian solution using the Lagrangian state 
vector is presented in the next chapter. 
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Figure 5.1: Local coordinate system in cell with node numbering and marker location 
at r = ( z,y,z ). 

5.1 Convection of the markers and integration of the 
source-terms 


At a given time-step in the Eulerian solution a set of markers is initialized at a chosen 
location in the Eulerian grid. Associated with the markers are arrays indicating in which 
cell each marker is located. The initialization of these arrays is provided by a brute force 
search through the whole domain as described in Appendix D. As the markers move 
through the Eulerian grid, these arrays are reset with the proper cell values by using a 
system of ‘neighbouring cells’ which will be described later. 

A cell-centered local coordinate system (£,»?, C) * s set U P ™ eac h ce ^ °f the compu- 
tational domain as shown in Figure 5.1. The eight tri-linear interpolating functions Ni 
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to N 8 are defined from cell comers to marker in the cell as 


*1 =1(1 -0(1 -^(1-0, 
^2 =1(1 + 0(1 -0(1-0, 
n 3 =1(1 + 0(1 + 0(1-0, 
n 4 =1(1 -0(1 + 0(1-0, 


N s =i(i-0( 1 - , ?)(i + C)> 
N e = |(i + 0(i -0(i + 0. 
Ni = 1(1 + 0(1 + 0(1 + 0 , 
■W* = s(i - 0(1 + */)(! + 0-- 


At each iteration, the local cell coordinates (£, 77 , C) f° r a marker axe determined by 
solving the implicit system 


8 

? = £-Mfe(f, V, 04 (5.2) 

*=1 

by a few Newton-Raphson iterations as described in Appendix E. r and f *. are the 
marker and the nodes positions, respectively. Usually, only three Newton-Raphson 
iterations are required. 


Any Eulerian function / (say the velocity or the entropy) defined at the cell nodes 
is transferred to the marker location in the cell by the tri-linear interpolation as 

8 

/( 0 » 7,0 = £^*( 0 * 7.0 fk- (5.3) 

fc=l 

For steady-state solutions, the Eulerian and the Lagrangian solutions do not need to be 
integrated in time using the same time- steps during the pseudo-unsteady convergence. 
A local time-step can be used for the markers integration since, for steady flows, the 
markers can convect at any speed along the stre amlin es. In unsteady flows, however, 
the Eulerian and Lagrangian time-steps have to be the same, since at each instant in 
time, the markers represent the current location and Lagrangian state vector of a given 
particle. In the following, the Eulerian and the Lagrangian time-steps will be denoted 
At and At/, respectively. 

When a time-step At/ is taken in the Lagrangian solution in order to advance the 
marker position and integrate the source- terms, the new marker position is found by 
using a predictor- corrector integration scheme as 

Tp = f\tl) + At/ 
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(5.4) 


f c = r{ti) + At/ v? +1 , 

r[ti -+- A t/) = — (r p I " c ^’ 

where f p and r c are the predicted and corrected marker positions, respectively. ££ is the 
Eulerian velocity at time t interpolated from the cell nodes to the old marker position 
f, and u? +1 is the Eulerian velocity at time t + At interpolated from the cell nodes to 
the predicted marker position r p . The interpolation from nodes to marker position is 
performed using Equation (5.3). In System (5.4), the Eulerian velocity fields at time t 
and t + At are used to find the new marker position. The predictor and the corrector 
step can also be based only on the Eulerian velocity field at time t. Both options are 
used depending on the Eulerian/Lagrangian interaction procedure as will be described 
in Section 5.2. 

The marker’s local coordinates have to be recalculated for the updated position 
of the marker f(ti + At/) by again solving System (5.2). If any of the f, tj, £ values 
falls out of the range (-1,1) but in the range (-2,2), the marker has moved out of 
the current cell but is located in one of the neighbouring cells. As mentioned before, 
associated with each cell is an array containing its 26 neighbouring cells so that when a 
marker moves out of a cell it can be readily relocated. Moreover, the locations of the 26 
neighbouring cells relative to the central cell are known, so that the values of £,r) and 
£ directly indicate in which of the 26 cells the marker is located. (In the case of the 
grid singularities mentioned in Section 3.3 and encountered in pipe flow calculations, 
the missing neighbouring cells are accounted for by assigning two pointers to the same 
cell). The marker local coordinates are now recomputed by again solving System (5.2) 
for the cell where the marker has been found. 

If the marker has moved more than one cell (any of the £, »7, C values are out of the 
range (-2,2)), a brute force search through the whole domain is required. This time- 
consuming procedure is avoided by specifying a Lagrangian time-step smellier than the 
time- step allowed by the CFL criterion in the Eulerian solution. The marker will then 
move over a distance of less them one cell in one Lagrangian time-step At/. However, in 
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steady flows with strong streamline curvature, as shown by the results from the swirling- 
flow cases, the CFL restriction does not allow for accurate streamline integrations. In 
these flows, the predicted marker position is used to scale the Lagrangian time-step such 
that the marker cannot cross the cell in less than two Lagrangian time-steps in any local 
cell coordinate direction. 


The domain boundaries are bordered by ‘wastebasket cells’, so that any marker 
exiting the domain is either eliminated or can be reinitialized somewhere else in the 
domain. In order to minimize the array sizes, a marker renumbering system is used to 
eliminate any exiting marker from the arrays and introduce new markers in the domain. 


In parallel with the trajectory integration described by System (5.4), the Lagrangian 
values of vorticity are computed at each marker by a predictor-corrector integration of 
the source-terms for the vorticity from Equation (4.3) as 


AuT p = At i (d? • V)v — d>(V * v ) — V x Vp 

Au T c — Ati (d? • V)tf - u5(V • v) - V x Vp 


n+l 


(5.5) 


for a compressible flow case. In turn, from Equation (4.6) for an incompressible flow 


Au? p = AU [(£} * V)fi]£, 

AuT c = AU [(d? - V)^ 1 . (5.6) 

Then the Lagrangian value of vorticity at the new marker location is found by 

${ti + Ati) = u>{ti) + i( Ad3p + AdJ c ). (5.7) 

In Systems (5.5) and (5.6), the quantities within brackets are determined at the old 
marker location t by using the Eulerian state vector at time t for the predictor step 
and at predicted marker location r p by using the Eulerian state vector at time t + At, 
As for the trajectory integration, depending on the Eulerian/Lagrangian interaction 
procedure, predicted and corrected values can also be dependent on the Eulerian state 
vector at time t only. 
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The determination of the vorticity and of the non-barotropic term V(1 / p) X Vp at 
the marker location (old or predicted) on the right-hand side of Systems (5.5) and (5.6) 
requires the definition of spatial derivatives. These are found directly for any Eulerian 
function / defined at the grid nodes by differentiating Equation (5.3) 

df 4-(dN k dj dNkdr, dN k 3C\ f f5 g) 

df " V d( dr drj df d( df ) k 

The metrics’ derivatives (d£/ df),{dr)l df),{dQI df) are directly found when using the 
Newton-Raphson procedure as described in Appendix E. The velocity derivatives com- 
puted by Equation (5.8) are required for the integration of the source-terms. In a 
compressible flow, the spatial derivatives of density and pressure have also to be com- 
puted. 

In addition to Equation (5.5), for a compressible flow, the entropy-related function 
5 is integrated along trajectories as 

5(tj + At,) = S(ti). (5.9) 


5.2 Eulerian/Lagrangian interaction 

In order to achieve an efficient interaction between the Eulerian and the Lagrangian 
solutions, two different strategies have been implemented for the different flow cases 
treated here. 


5.2.1 Downstream integration of trajectories 

The first and more straight forward option consists in positioning the markers at the 
inlet of the computational domain and tracing them as they convect downstream with 
the local Eulerian flow. At each time-step along their respective trajectories, the new 
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Lagrangian state vectors are obtained by integration of the source-terms and the markers 
are used to locally correct the Eulerian flow. The initial values for the Lagrangian state 
vectors are provided by linear interpolation from the Eulerian flow at the inlet of the 
computational domain where the inlet boundary conditions are known. At this location, 
the Eulerian solution does not significantly suffer from numerical diffusion due to the 
proximity of the boundary conditions. For example, the vorticity field at the inlet of a 
pipe will be function of the specified inlet velocity profile. 

Figure 5.2b) shows a vortex connecting (perpendicular to its axis of rotation) through 
a contraction. For this unsteady calculation, the markers are initially positioned in the 
vortical region of the flow (near the vortex core) and then convect downstream together 
with the vortex. In this unsteady test case, the time-step A ti used for the integration 
of the marker trajectories must be the identical to the Eulerian At. Figure 5.2a) shows 
the sequence of interactions between the Eulerian and the Lagrangian solutions. First 
the Eulerian solution at time t and the Eulerian solution at time t -f A t are used to 
compute the new marker position following Equations (5.4). The new values of vorticity 
and entropy are computed using Equations (5.7) and (5.9) for a compressible flow and 
Equation (5.7) for an incompressible flow. The new Lagrangian state vector is then used 
to correct the Eulerian state vector for each cell containing a marker at time t -f At. 

Figure 5.2c) represents a schematic for a steady swirling flow through a pipe. The 
markers are injected from the inlet of the domain at regular time intervals in such a way 
as to have approximately one marker per cell in the vortical regions of the domain. As 
markers are converted downstream and exit the domain, it is necessary to reinject them 
at the inlet. The trajectories and the source-term integration are computed in the same 
m a nn er as for the unsteady process but now since the Eulerian solver uses a pseudo- 
time marching, the Lagrangian time-step does not need to equal the Eulerian time-step. 
The Lagrangian time-step is now fixed by accuracy considerations in the trajectory 
integration. The interaction between the Eulerian and the Lagrangian solvers is again 
described by Figure 5.2a). 
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The downstream integration of the trajectories presents the advantage that no a 
priori knowledge of the location or structure of the flow features is necessary. For an 
unsteady flow, the markers axe simply positioned in the features area at the initial time. 
The use of identical Lagrangian and Eulerian time-steps ensure the presence of markers 
in the vortex area for each time. For the steady swirling flow in a pipe, the markers can 
be positioned in the entire inlet cross-section. In this case, however, it can be inferred 
that placing markers only near the wall region of the pipe at the inlet is sufficient to 
trace the strong vorticity regions downstream. 

On the other hand, for this last flow example, the markers are subjected to a strong 
redistribution during their convection, such that an even correction of the Eulerian 
solution can not be ensured. In order get a more even distribution of markers in the 
areas where a correction is required, a second Eulerian/Lagrangian interaction strategy 
is used. 


5.2.2 Upstream integration of trajectories 

The second interaction strategy is illustrated in Figure 5.2d) and (5.2e). For the same 
steady swirling flow through a pipe, the markers are initially placed at the center of each 
cell and the trajectories and source-terms are integrated upstream until the markers 
reach the inlet. The time-step for the integration is again fixed for each marker by 
accuracy considerations in the trajectory integration. At the inlet, the Lagrangian state 
vectors are found by linearly interpolating the inlet (non-diflused) Eulerian solution to 
the markers location. By adding the integrated source-terms to the inlet state vectors 
the marker state vector at the center of each cell is determined. The correction procedure 
can then take place evenly from the center of each cell. Figure 5. 2d) describes the 
interaction between the Lagrangian and the Eulerian solvers. The integration of the 
source-terms requires interpolations from the Eulerian state vector at time t. Once the 
value at the center of each cell is determined, the correction of the Eulerian solution is 
performed. Then the Lax-WendrofF algorithm is applied on the corrected solution. At 
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this point, the trajectories could be computed again, but it has been found that iterating 
the correction step and the Lax- WendrofF algorithm, while keeping the Lagrangian state 
vector and the marker trajectories fixed, leads to a faster convergence. On the other 
hand, to ensure stability, the trajectories should be recomputed before the Eulerian 
solution changes are too large. Also, the trajectories do not have to be traced backwards 
to the inlet, but can be stopped at any previous upstream location. If the trajectories 
are stopped before the inlet, the corrections will be smaller and the convergence will take 
more time since the full correction is obtained only when the cells where the markers 
stop have been fully corrected. Again, in order to minimi ze CPU time, the markers can 
be placed only in the flow regions where corrections are required. 

For an unsteady flow solution using the upstream integration of the trajectories 
means that the trajectories have to be recomputed at each time-step. Again, the tra- 
jectories can be integrated backwards until they reach the inlet or be integrated over 
only a few time-steps. At any backward position, though, the Eulerian solution at that 
particular time is required to compute the source-terms and the trajectories. Therefore, 
in order to minimize the number of Eulerian solutions to store, the trajectories can 
be integrated backwards over only one time-step. For unsteady flows, the downstream 
integration of trajectories is, therefore, simpler and much more efficient. 

In summary, the downstream integration is the more straight forward method, es- 
pecially for unsteady flows. For steady flows, the use of the upstream integration of 
the trajectories allows a spatially-even correction of the Eulerian solution. This has 
met with more success for the calculation of swirling flows and secondary flows. Both 
trajectory integrations are further discussed in relation with the markers positioning in 
the flow. 
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d) 


Figure 5.2: Eulerian/Lagrangian interaction procedures for downstream integration of 
the trajectories: a) schematic of Eulerian/Lagrangian interaction, b ) unsteady vortex 
convection in contraction, c) steady swirling flow in pipe, for upstream integration of 
the trajectories: d ) schematic of Eulerian/Lagrangian interaction, e) steady swirling 
flow in a pipe. 
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5.3 Positioning of the markers in the flow 


Since each marker corrects the Eulerian solution locally, the markers need to be located 
only in the regions of the flow where a correction is required. This property is particu- 
larly useful in terms of CPU time reduction when the flow features are concentrated in 
an otherwise smooth background flowfield. In contrast, the ‘cloud- in- cell’ technique of 
Baker [5] uses finite-size vortices which influence the entire flowfield and whose geometry 
and circulation must be known. 

When using a downstream integration, the built-in convection properties of the 
Lagrangian technique allow for the steady or unsteady tracing of flow features without 
an a priori knowledge of their placement or structure. In the upstream integration of 
the trajectories, the area covered by the flow features must be broadly known before. 
It is usually possible to determine the general area where corrections are required from 
a basic knowledge of the flow. However, in some cases like the secondary flow in a pipe, 
the features of interest are so dispersed in the flow that markers have to cover almost 
the entire flow area. 

If the markers are convecting downstream and locally correcting the Eulerian flow, 
weighting factors are required when distributing corrections to the grid nodes because 
at any time-step the markers are not distributed uniformly with respect to the Eulerian 
grid, (each node may not be influenced by the same number of markers and the distances 
from the nodes to the markers are different). Moreover the divergence /convergence of 
the trajectories can create holes in the markers distribution. In the case illustrated in 
Figure 5.3a), because of the divergence of the stre amlin es, the markers are redistributed 
as they convect from the inlet downstream. This results in some nodes in the field 
not being influenced by any markers, whereas other nodes are influenced by several 
markers. When the correction is large enough, this uneven correction applied to the 
Eulerian scheme can lead to numerical instability. 

The use of the second interaction option where the markers are placed at the center 
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of each cell and convect upstream eliminates the need for weighting functions. Since the 
distribution of the corrections is done from the cell centers, each grid node is influenced 
by the same number of markers disposed uniformly. The problem of redistribution of 
the markers during their convection is also eliminated as presented in Figure 5.3b). 
Nevertheless, as pointed out before, the choice of upstream integration is less suitable 
for unsteady flow calculations. 

The assumption that the Lagrangian state vector is piecewise constant in each cell 
makes the representation of strong gradients inaccurate on very coarse grids. This case is 
illustrated in Figure 5.3c) and 5.3d) for the downstream and upstream integration of the 
trajectories, respectively. In Figure 5.3c) where the markers are convecting downstream, 
using the weighting factors results in an average correction for the field cells. This is 
clearly inaccurate if the gradient between the Lagrangian state vectors is high. In 
the case of upstream integration represented in Figure 5.3d), the presence of only one 
marker in each field cell is ensured. However, because of the lack of grid resolution 
and the presence of only one marker per cell in the field region, when the markers 
trajectories are diverging (in an upstream integration), the final distribution of markers 
will be sparse in the inlet region as shown in Figure 5.3d). The inlet information 
between the two stre amlin es is not ‘seen’ or transported by the markers. Both problems 
illustrated in Figures 5.3c) and 5.3d) are linked to a lack of grid resolution. The 
problem of high vorticity gradients does not appear in the Eulerian solution alone, since 
numerical diffusion results in lower vorticity gradients in the field. The Lagrangian 
solution, however, is immun e to numerical diffusion and presents stronger vorticity 
gradients. 

The different flow cases treated here reflect the need for distinct marker initial lo- 
cations and trajectory integrations. For the unsteady vortex convection, few markers 
are required in the flow and the markers are attached to the flow feature as it convects 
downstream. For the steady swirling flow in a pipe, many more markers are required 
and are initialized near the inlet cross-section. However, since a strong redistribution 
occurs in a downstream integration, one can also attempt an upstream integration of 
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Inlet Field 



c ) 


d) 


Figure 5.3: Markers redistribution during downstream or upstream convection leads to 
a) lack of correction of the Eulerian solution, b) even correction, c) average correction 
in cell, and d) sparse distribution of markers in the inlet region leading to inaccurate 
representation of inlet flow values by the Lagrangian markers. 
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the markers by initially locating a marker at the center of each cell. This last procedure 
ensures an even correction distribution, but the downstream integration of the trajec- 
tories remains much simpler for unsteady flow calculations. The lack of grid resolution 
can create problems when tracking the Lagrangian vorticity upstream or downstream. 
Section 9.2 addresses this problem by using a pseudo-viscous term in the Helmholtz 
equation in order to smooth out the strong vorticity gradients of the flow. 
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Chapter 6 


Correction Procedure 


This section describes the corrections by which the Lagrangian solution induces 
changes in the Eulerian state vector. Since the vorticity and entropy are accurately 
described by the Lagrangian equations, these quantities are used to correct locally the 
Eulerian solution. For an incompressible flow the vorticity at the marker location is 
used to locally correct the Eulerian solution. An entropy correction is used in the 
compressible flow cases in addition to the vorticity correction. 

The main difference in terms of correction when using the downstream or upstream 
trajectory integrations described in the previous section is the location of the markers 
in the field. With a downstream integration the markers are not placed unif ormly in 
the flow. Therefore, the use of weighting factors is required when interpolating the 
flowfield values from the grid nodes or from the cell centers to the location of the 
marker and vice-versa. When using an upstream trajectory integration the weighting 
functions are not required in the correction procedure since each correction occurs at 
the initial location of the markers, i.e. from the center of each cell where a marker has 
been placed. The correction procedure is described here after for a non- unif orm marker 
distribution (or a downstream trajectory integration). The simplifications for a uniform 
marker distribution (upstream trajectory integration) are also mentioned. 
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6.1 Entropy correction 


For a non-uniform marker distribution in the field, the Eulerian value of the entropy- 
related function S c defined by Equation (4.4) is found at the marker location by linear 
interpolation from the nodes according to Equation (5.3) as 

8 

S e = Y,NkS k . (6.1) 

k—1 

This value is then compared to the Lagrangian value of the entropy-related function 5/ 
at the marker location. The difference 


AS = Si - S c (6.2) 

is considered as an error in the Eulerian solution at the location of the marker. This 
error is then distributed back to the nodes as illustrated in Figure 6.1. When using a 
downstream trajectory integration, each node may be influenced by a different number 

of markers, so that the resulting change in 5 at a node is found by a weighted average 

over all the markers influencing that node (i.e. all the markers located in cells to which 
that node belongs) 

E-fttC - * (6 ' 3) 

If the markers are located at the cell centers, Equations (6.1) and (6.3) simply become 

Se = l £ S ( 6 - 4 ) 

0 k~l 

and 

ss = l £ A5 m , (6.5) 

” m= 1 

respectively. The change in the 5 value at a node is then translated into changes in the 
Eulerian flow field variables. 


It is necessary to fix four of the five flow variables in order to induce a change in 
entropy in the Eulerian solution. If the conservation variables ptx, pv, pw and the pressure 


88 




Figure 6.1: Distribution of the error in entropy-related function AS from markers to 
nodes using weighted average. 

are kept constant on the grounds that these elliptic quantities are well represented on 
the Eulerian grid, then the change in S is translated into changes in the conservative 
variables as 


6p = 

-Ass 

7 S 

(6.6) 

S(pu) = 

0 

(6.7) 

6(pv) - 

0 

(6.8) 

S(pw) = 

0 

(6.9) 

6{pe 0 ) = 

u 2 +v 2 +w 2 

2 ‘f 

(6.10) 


Another option, si m i l ar to the approach used by Giles for his entropy smoothing 
in [23], consists in keeping the vorticity constant during the entropy correction since 
the correction in vorticity is the object of a separate treatment. This is achieved by 


89 


updating the density and energy and keeping the velocity and pressure constant. Then 


6p = 

--L6S 

(6.11) 

6{pu) = 

u8p 

(6.12) 

8{pv) = 

v8p 

(6.13) 

8(pw) = 

w8p 

(6.14) 

6(pe 0 ) = 

u 2 + v 2 + w 2 

2 Sf 

(6.15) 


In practice, both distribution methods give the same results. 


6.2 Vorticity correction 


6.2.1 Vorticity error at cell centers 


The Eulerian vorticity 3 is first computed at the center of each cell by assuming a 
cellwise constant vorticity. Then 


/ 3dV = 3V c 
Jv 


(6.16) 


where the integral is performed over the cell volume V c . Using Stokes theorem the 
vorticity is 

3 = f 3dV = tt / (V X v)dV = f (n x v)dS (6.17) 

V c Jv V c Jv K Js 

where the surface integral is performed over the cell faces and n is the surface unit 
normal vector. Thus, s ummin g over the six cell faces, the components of vorticity in 
Cartesian coordinates are 

6 




— — Sy—V S Z )i 


i-l 
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( 6 . 18 ) 


1 

u y = TF Sz - ™ S *)” 

Vc i = 1 

1 6 

= 7F - u S y)i 

c »=i 

On each face, the values for u,t? and w are obtained from averages over the velocities 
at the four corner nodes. 

In the case of a non-uniform marker distribution, the Eulerian vorticity is then 
interpolated to the location of the marker by using a second local coordinate system 
based on the eight cell centers that the marker is the nearest to (these 8 
cells are easily determined by using the array containing the 26 cell neighbours and the 
marker local coordinates (f , 77, ()). The £, 77, ( directly define the (£*, 77*, (*) coordinates 
as shown by the example of Figure 6.2 (represented in a 2-D case for clarity). 

The 8 tri-linear interpolating functions 77% (*) to 77*, C) are used to 

interpolate the Eulerian vorticity from the eight cell centers to the marker location as 

8 

= ( 6 . 19 ) 

C= 1 

This value is then compared to the Lagrangian value of vorticity at the marker location 
( 3 ;. The difference 


Au; = u;/ - Q e (6.20) 

is taken as an error in the Eulerian solution at the location of the marker. This error is 
first distributed back to the center of the eight cells by using a weighted average over 
ail the markers influencing this cell center. 

(6 - 21) 

All the markers influencing a cell center are contained in a volume based on the eight 
neighbouring cell centers as illustrated by the dashed area in Figure 6 . 2 . If the markers 
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Figure 6.2: Local coordinate system (£*,??*, C) centered on nodes (represented for the 
solid line marker) and vorticity error distribution from markers to cell centers. 

are located at the cell centers, Equations (6.19) and (6.21) become 

= u? c , (6.22) 

and 

8Q C — Au>, (6.23) 

respectively. Once the change in vorticity at a cell center is determined, it is translated 
into changes in the Eulerian velocity vectors at the comer nodes of the cell as described 
in the next section. 
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6.2.2 Distribution of vorticity error 


The error in vorticity in a cell is translated into changes in the velocity components 
at the cell nodes. In a 2-D problem, the error in vorticity implies a circulation error 
around the cell perimeter. The correction is then distributed evenly between the eight 
velocity components at the four cell nodes as described in [20]. In three dimensions, the 
vorticity error is a vector, but a similar procedure can be applied for each of the six cell 
surfaces. 

The error in vorticity Au; at the cell center is used to compute an error in circulation 
on each face of the cell as 


AT = 5, (6.24) 

where S is the surface vector of the face. Figure 6.3 shows the errors in circulation for the 
faces 1 4 and 6 (the nomenclature for the face numbering is described in Appendix B). 
On each face, the error in circulation is then used to find velocity corrections at the 
four comer nodes as shown in Figure 6.3a). A local coordinate system (< 7 , r) and a local 
node numbering from 1/ to 4/ is defined on each face as sketched in Figure 6.3b) for 
face number 1. The location vector f = ( x , t/, z) is used to define the unit vectors along 
the local coordinates as 

a= r 2f - r lf + r 3f ~ r 4f f= r 4/ - f lf + f 3/ - f 2/ 

11^2/ - ?\f + r zj - r Af || 5 ||r 4/ - r lf + f zf - r 2 f\\' 

The local coordinates a and r are then written as 

a ~ a • (f - f c ), r = f • (r - f c ), (6.26) 

where f c refers to the location vector of the geometric center of the face. 


The error in circulation can be written as an integral over the face perimeter as 

AT = A udx + A vdy + Awdz = Au a da + A u T dr, (6.27) 

where A u, Av, Aw are corrections to the u, v and w components of velocity, respectively. 
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Figure 6.3: a) Velocity correction (contribution from faces 1 4 and 6 only), and b) local 
coordinate system < 7 , t and local node numbering for face number 1. 

Auj, , Atir are corrections to the u a ,u T components of velocity expressed in the local 
coordinate system. Using trapezoidal integration, the error in circulation is written as 

2 Ar = + Au„ 2/)(^2/ - °1/) + (Attrl/ + AUt2/)(t2/ — Tif) + 

(Au < t2 / + Au„3/)(ff3/ - 02 /) + (AtiT-2 / + Aii T 3/)(r 3 / - r 2 /) + 

(A U aZf + Au„4/)((T4/ - cr 3 f) + (AUt 3 f + AUt4/)(T4/ ~ T 3 /) + 

(Au„ 4 / + Au <r i/)(o'i/ - cr 4 /) + (A V~rAj + Ati T i/)(r 1 / - r 4 /) 

= (Att„i/ - A U ff 3 /)(<T 2 / - 04/) + (Au a2 f - Au a 4 f)(<T 3 f - 01/) + 

(Au t i/ - Au t3 /)(t 2 / - T 4 f) + (A Ur2f - Au t4 /)(t 3/ - T lf ). (6.28) 

By choice, the error in circulation is translated into homogeneous velocity corrections 
for each of the eight velocity components at the four comer nodes. Then 

AtXoyr = Au„i/ = A U a3 f = -At l a3 f — -Au ff4 / 
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— AUrlf — A U t2 / = A U T zf — -AtX T 4 /. 


(6.29) 


Thus, the correction to the components u a ,u T is found as 


A V'cr^r — 


AT 


where P — <73 y — <Ji/ + < 73 / — a^f + r 4 ^ — 73 / + r$f — t^. (6.30) 


The transformation from the local face coordinates to the Cartesian coordinates in 
Equation (6.27) gives the velocity corrections A u, Au, Ait; 


Ar = (j) Au c da + A u T dr 

= /a + 3 ^ + £f*) + / 4 M§J< 1 * + §J<*»+ §J*) 

= /(*«,§; + A«,g)<fe + /(A»„|i + + /(A U „^ + 

= £ Au dz + At? dy + An; dz, 

so that 


where 


An 

An 

An; 


A dcr K dr 

An^— + An r — — , 
dz oz 

dcr dr 

Ati^— + Au r — , 

A dcr dr 

Au "«I + A ^«? 


( 6 . 31 ) 


da 

x 2f ~ x lf + x 3 / ~ x 4 f 

dr 

x 4f ~ X 1 / + x 3f ~ x 2f 

fa 

11*2/ - *1/ + *3 j - 7\t/|| ’ 


11*4/ - *1/ + *3/ — 7=2/11 ’ 

da 

V2f - yif + y3f - 2 / 4 / 

dr 

2/4/ - 2/1/ + 2/3/ ~ 2/2 / 

dy 

11*2/ - *1 f + *3/ - 7=4/11 ’ 

dy 

11*4/ -ri/ + f 3/ -fj/ll’ 

da 

z 2 f ~ z lf + z 3f ~ z 4 / 

dr 

z 4 / — z lf + z 3 f ~ z 2f 

dz 

11^*2/ - *i/ + r 3/ - r 4/ ||’ 

J~z ~ 

11*4/ ~ *1/ + *3/ “ *2/lf 


At a node, the contributions to the velocity change from each face are added to give the 
final velocity correction since, if the faces are perpendicular to each other, the velocity 
correction on one face does not affect the circulation on another face. Each node receives 
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Figure 6.4: Solid-body rotation distribution of error in vorticity to velocity components. 

contributions from three faces for a given cell. The final node velocity correction is then 
found by adding the contributions from the eight cells surrounding the node. 

This velocity correction method is analogous to the use of a solid-body rotation 
assumption where the change in vorticity at a cell center induces changes in velocity 
components at the 8 nodes of the cell as 

6v = 6Sl c X f c , (6.32) 

where the change in angular speed 6£l c in the cell is defined as half the vorticity error 
in the cell 

*n c = (6.33) 

and f c is the vector connecting the node to the center of the cell as illustrated in 
Figure 6.4. The analogy between the previous distribution procedure and the solid- 
body correction does not hold if the cell faces are of different shapes. Using the pre- 
viously described correction distribution ensures identical corrections for faces of same 
(area/perimeter) ratio, independently of the face shape as illustrated in Figure 6.5. 




Figure 6.5: Distribution procedure ensures identical velocity correction for faces of 
same (area/perimeter) ratio, independently of face shape as opposed to the solid-body 
rotation distribution (dotted line vectors). 

Instead, the correction obtained through a solid-body correction does not lead to an 
identical correction in every velocity component as indicated by the dotted arrows. In 
the case of a sheared cell, since the local coordinate system is not orthogonal, the correc- 
tions are different from the orthogonal case. However, the right correction in circulation 
and vorticity are still enforced. 


For an incompressible flow, summing the contributions from each face described by 
Equation (6.31) and averaging over the eight surrounding cells, together with 

= 0, (6.34) 

give the changes in the Eulerian conservative variables due to the vorticity correction. 
For a compressible flow, the resulting changes in the conservative Eulerian variables are 
found by 


Sp = 

o, 

(6.35) 

6{pu) = 

pSu, 

(6.36) 

6(pv) = 

pSv, 

(6.37) 

S(pw) = 

p Sw, 

(6.38) 


97 





6(pe 0 ) = 0. 


(6.39) 


The actual implementation of the correction procedure is performed in two steps. First, 
the contribution from each surrounding face to a node is computed for a unit value 
of AT on the face (using Equations (6.30) and (6.31) with AT = 1). For example, 
in Figure 6.3 node 2 receives contributions from faces 1, 4 and 6. By using a unit 
value for AT, this step reduces to a simple metrics calculation required only once at the 
beginning of the computation. Then, at each iteration in time, the error in vorticity 
at the cell center is determined by Equation (6.21) or (6.23) and the the contributions 
for a unit circulation are multiplied by the actual value of AT for each face. Then the 
contributions of the three faces are summed for each node. The velocity correction at 
a node is then found from by adding the corrections due to the eight surrounding cells. 
Since this last step involves a loop over the cells, it is fully vectorizable. 


6.3 Boundary conditions for velocity correction 


The nodes located on the limits of the domain receive contributions from less than eight 
cells and, therefore, require special treatment. The imposition of boundary conditions 
is performed during the first step process for nodes located on do main boundaries and 
results in modifications of the contributions to these particular nodes. 

In this work, four types of domain boundaries are present : inlet, exit, wall and 
symmetry-plane. The correction of velocity for the nodes located on an inlet boundary 
is set to zero, since the Eulerian values are exactly set through the inlet boundary 
condition. The nodes located on the exit surface miss the contributions from the cells 
placed downstre am of the exit surface. Pseudo-cells are defined with the same error 
in vorticity at their center than the one for cells lying inside the domain as illustrated 
in Figure 6.6a). This is the method used for the flow over a three-dimensional wing 
where the trailing vortex axis crosses the exit boundary almost perpendicularly. For 
the secondary flow calculation in a circular pipe, the axis of the secondary vortex crosses 
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Figure 6.6: Correction in vorticity for pseudo-cells placed at a) an exit surface, b) a wall 
surface, and c) a symmetry surface. 

the exit boundary with a smaller angle and a better approach is to use the same type of 
exit boundary condition as the Eulerian solver. The 1-D characteristics theory is then 
applied to a vector formed by the corrections in velocity components. 


The imposition of wall and symmetry-plane boundary conditions uses pseudo- cells 
as sketched in Figures 6.6b) and 6.6c), respectively. 

If C is the local coordinate normal to the exit, wall or symmetry-plane, then 


&&cxit 

= (Aw ( ,Aw„Aw ( ), 

(6.40) 



(6.41) 


= (-Au> { , -Au n , Au> ( ). 

(6.42) 


The corresponding velocity corrections are sketched in Figure 6.7 for faces 6 , 4 and 
1, respectively. The additional velocity corrections due to the presence of the pseudo- 
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cells are represented as dashed vectors. Figure 6.7a) shows the case of an exit boundary 
condition whereas the case of a wall or symmetry-plane boundary condition is presented 
in Figure 6.7b). 


6.4 Discretisation error for velocity correction 


The discretisation error obtained by using the Lagrangian correction on a finite-size 
grid is estimated here for the case of a linear vorticity profile (say a Poiseuille flow in a 
channel). The vorticity is written as a linear function of channel height z as 

u > = kz, (6.43) 

where k is a proportionality factor. In the Lagrangian technique, the vorticity is assumed 
to be a constant within each cell and takes the value at the cell center. By expanding 
the vorticity as a function of the channel height into a Taylor expansion to the first 
order, the discretisation error e u in a cell is 

e u = tt-Az ~ k l, (6.44) 

az 

where Az is the distance to the cell center and l is a representative cell length. Using 
Equation (6.30), the velocity correction An is related to the vorticity correction as 

Au = — = ~ Aw Z. (6.45) 

P P 

Thus, the discretisation error in velocity correction e u is 

e u ~ e u l (6.46) 

and using Equation (6.44), the error in velocity is shown to be second-order on the 
finite-size grid 

e u ~ k l 2 . (6.47) 
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b) 



Figure 6.7: Contribution to velocity corrections from faces 6, 4 and 1 in the case of a) 
exit boundary condition, and b) wall or symmetry-plane boundary condition (additional 
contributions from pseudo-cells tire represented as dashed line vectors). 
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6.5 Iterative procedure for the correction of vorticity 


The previous section’s subject was the translation of an error in vorticity, defined at the 
cell center, into errors in circulation around each face of the cell, and finally into velocity 
corrections at the comer nodes. This section shows why a ‘recursive’ correction process 
is required in order to obtain the correct circulation around each cell. The different 
steps of the iterative correction procedure are presented as well as a numerical study of 
the process’s convergence rate for a simple two-dimensional case. 

As will be shown, the vorticity correction requires an iterative procedure. For steady 
flows, the vorticity correction is converging along with the pseudo time-marching Eu- 
lerian scheme and only one step is usually required for the vorticity correction at each 
Eulerian iteration. In the case of an unsteady flow, the vorticity correction should be 
iterated until convergence for each Eulerian time-step. However, using only one step of 
the vorticity correction at each Eulerian iteration has given satisfactory results. 

In this study, the correction is applied recursively on an initial solution without 
taking any Eulerian steps. However, the accuracy of the Eulerian solution enters into 
the analysis of the convergence rate of the combined Eulerian/Lagrangian scheme. 

Figure 6.8 shows a two-dimensional numerical domain with three cells in each direc- 
tion. The initial solution is such that the circulation in each cell is zero. The desired 
circulation in each cell is assumed to be unity. Therefore, a uniform circulation cor- 
rection is imposed for each cell (Ar = 1). In the first step, velocity corrections are 
determined for each cell such as to satisfy the circulation correction according to Sec- 
tion 6.2.2. The resulting velocity correction at a node is then found by adding the 
contributions from the surrounding cells. When all the contributions are added, the 
circulation in each cell is modified. Hence, a recursive correction is required. In the 
second step, the modified circulation in each cell is compared to the desired circulation 
in order to define a new circulation error. This is in turn translated into new velocity 
corrections at the comer nodes. The procedure is iterated until the correct circulation 
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in each cell is obtained. 


Figure 6.9a) shows the maximum error in circulation Ar max in the domain as a 
function of the iteration number n for six grid sizes (2 x 2 to 33 X 33 nodes), each grid 
having twice as many cells in each direction than the previous one. ATm ax is normalized 
by the initial constant error in the domain Ar initia ;. The number of iterations that 
the information takes to travel from the boundaries to the center cell is shown by a 
constant error in circulation during the first iterations (the number of iterations with 
constant error corresponds to half the domain size in one computational direction). The 
convergence rate is slower with increased grid size. However, each time the grid size is 
increased, the initial error in the domain is reduced by a factor 4 since the Lax- Wendroff 
scheme is second-order accurate. 

In Figure 6.9b), the ratio of the maximum error in circulation between two con- 
secutive iterations is shown to converge to a value dependent only on the grid size. 
Figure 6.10 displays the behavior of the difference between these values and unity (1 
represents no convergence) as a function of the mesh size h, plotted in a logarithmic 
scale. This measure of ‘how slowly the correction procedure is converging’, is shown 
to vary as a second- order function of the mesh size. However, the initial error given 
by the Eulerian scheme is also a second-order function of the grid size. Thus, as the 
convergence rate of the iterative procedure goes to zero at the limit of an infinite grid 
resolution, the truncation errors due to the Lax- Wendroff scheme vanish too. Thus, the 
resulting solution is still consistent. 

The extension to three dimensions is straight forward. As mentioned earlier, the 
velocity corrections on one face do not induce any circulation on another face if the 
faces are perpendicular to each other. Thus, the results for a two-dimensional case 
remain unchanged in three dimensions. (However, the velocity corrections have to be 
divided by a factor 2 for the ‘interior 5 faces since when perfor min g the three-dimensional 
vorticity correction in a unstructured fashion by a loop over the cells, each face (except 
the boundary faces) is taken twice into account). Also, the convergence rate remains 
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unchanged if the cells faces are rectangular instead of square. 


In the steady-state calculations, at each iteration the corrections can be multiplied by 
a factor sm all er than 1 in order to limit the perturbations to the Eulerian solution. In the 
flow over a 3-D wing, the magnitude of the correction in vorticity for the trailing vortex 
is of the same order as the basic vorticity in the flow. If the corrections are limited, the 
convergence rate to the correct circulation around each cell is slower and more iterations 
are required between the Eulerian and the Lagrangian solutions. Another option consists 
of under-rela xin g the circulation corrections by using the ‘old’ circulation correction 

Ar n+1 = AF* + Rf{AT n+1 - Ar n ), (6.48) 

where Rf is a relaxation factor smaller than 1 and n and n + 1 refer to the old and 
present iterations of the correction procedure, respectively. 
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Figure 6.9: Convergence of recursive vorticity correction: a) ma x imu m error in circula- 
tion in domain (referenced to initial error in circulation) as a function of the iteration 
number n and b) ratio of the maximum circulation error between two consecutive iter- 
ations. 
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Figure 6.10: Measure of ‘how slowly the iterative procedure is converging’ as a function 
of the grid size. 




Chapter 7 


Vortex Preservation Test Case 


The convection of a Lamb vortex in a three-dimensional uniform flow is computed 
as a vortex preservation test case in a compressible unsteady flow field. The vortex 
convects perpendicularly to its axis of rotation along a straight channel over a distance 
of approximately 50 core radii. The tangential velocity in the frame of reference of 
the Lamb vortex is given by 

f 7 ' 1 ) 

where T is the vortex circulation, a is the vortex core radius and r is the distance from 
the center of the vortex. The initial Eulerian and Lagrangian state vectors are obtained 
by superimposing the vortex on the uniform background flow. The entropy is found as 
a function of the radius by numerically integrating the radial entropy gradient given by 
Crocco’s relation 

(*- - 5*0 ™ - 'IT ^ <»> 

where v$ is given by Equation (7.1) and the stagnation enthalpy h is assumed con- 
stant in the vortex frame of reference. Once the entropy-related function S defined by 
Equation (4.4) is found, the density and pressure follow from 


* = 

(7.3) 

7 - 1 /. 1 2x 

P = P ( h °* ~ 2 V «) 

(7.4) 
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Figure 7.1: Computational grid for vortex preservation case (129 X 17 X 9 nodes). 

A unif orm flow Mach number of 0.7 and a value for T/a of twice the freestream stag- 
nation speed of sound are used in this calculation. The resulting maximum tangential 
Mach number in the frame of reference of the vortex is approximately 0.2. 

The grid used for the calculation is composed of 129 X 17 X 9 nodes and is shown 
in Figure 7.1. The channel width and height are 1/8 and 1/16 of the channel length, 
respectively. The mesh spacing is uniform and equal in each direction and chosen such 
as to have four cells across the vortex core diameter, hence 

*1 = *1 = = 2. (7.5) 

a a a 

A ‘no-flux’ boundary condition is imposed on the top and bottom walls of the channel. 
On the inlet, exit and side surfaces of the channel, the exact solution of the convecting 
vortex is imposed. This allows a minimi zation of the size of the computational domain 
(the half-width of the channel measures only four vortex core radii). 

The results are presented on a mesh surface at channel mid-height (the solution 
is independent of the channel height). The initial and final marker distributions for 
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Figure 7.2: a) Initial distribution of markers with core size and b) final distribution. 

this calculation are represented in Figures 7.2a) and 7.2b), respectively. The markers 
are initially placed in the vortex neighbourhood with one marker per cell. Additional 
markers are placed on the vortex axis (only 17 markers are located in the vortex core 
in a 2-D plane). The total number of markers for this calculation is 1160. The Eulerian 
solution and the Eulerian/Lagrangian solution are performed on the same grid. For this 
unsteady flow case, the same time-steps are used in the Eulerian and Lagrangian integra- 
tion scheme, hence the markers convect at the same speed as the vortex. The Lagrangian 
state vector consists of the position, the vorticity and the entropy-related function 5. 
The trajectories are integrated downstream and the computation of the source-terms 
for the vorticity gives no contribution due to the absence of three-dimensional effects. 
The interaction sequence between the Eulerian and the Lagran gian solvers follows the 
procedure described in Section 5.2.1. Only one correction step is performed at each 
Eulerian iteration at the current location of the markers. 

Figures 7.3a), 7.4a) and 7.5a) show the initial pressure, vorticity and 5 function 
contours, respectively, when the vortex is located at the inlet of the channel (all flow 
values are referenced as described in Section 2.3). After the vortex has traveled approxi- 
mately 53 core radii, the fined pressure, vorticity and 5 function contours are represented 
in Figures 7.3b), 7.4b) and 7.5b) when using the Eulerian formulation alone and in 
Figures 7.3c), 7.4c) and 7.5c) when using the Eulerian/Lagrangian scheme. The initial 
velocity vectors (after subtraction of the convection speed) are plotted in Figure 7.6a). 
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The vortex at the end of the channel is represented in Figures 7.6b) and 7.6c) for the 
Eulerian solution and the Eulerian/Lagrangian solution, respectively. 

A strong diffusion is observed in the Eulerian solution, whereas the addition of the 
Lagrangian correction is successful in preserving the vortex structure and intensity. 

This is substantiated by looking at the velocity (minus convection speed) and pres- 
sure profiles across the vortex (after traveling 53 core radii) represented in Figures 7.7 
and 7.8. Whereas the diffusion of the Eulerian solution is obvious, the vortex core 
radius and vortex strength remain almost unchanged when using the Lagrangian cor- 
rection procedure. 


Figure 7.9 shows a measure of the numerical diffusion as the decay rate of the pressure 
coefficient Cp v at the center of the vortex 


Cp v — 


Pmin Poo 



(7.6) 


where Pmin is the pressure at the center of the vortex and poo is the pressure of the 
background potential flow. The exact solution corresponds to a constant Cp v coefficient. 
For the same grid density, the standard Eulerian scheme leads to a strong diffusion 
of the vortex during its convection whereas a substantial improvement in the vortex 
preservation is obtained when using the combined Eulerian/Lagrangian scheme. The 
ragged aspect of the Cp v curve (computed as a cell-average value) when using the 
Lagrangian procedure is due to the vortex passing over the Eulerian grid. The difference 
between the Eulerian /Lagrangian curve and the exact curve is due to discretization 
errors when representing the vorticity and entropy of the vortex by a finite amount of 
markers placed at the cell centers. 


If the accuracy of the solution is measured by monitoring the pressure losses at the 
center of the vortex, the Eulerian/Lagrangian scheme leads to a solution approximately 
4 times more accurate than the Eulerian solution alone after the vortex has traveled 53 
core radii. In order to gain a factor 4 in accuracy a second-order accurate scheme such 
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as the Ni-Lax-WendrofF scheme would require a grid twice as fine in every direction. 
This would lead to an increase in CPU of a factor 16 (a factor 2 for every spatial 
direction times a factor 2 for the corresponding decrease in time-step). In comparison, 
the present Eulerian/Lagrangian scheme requires only ~ 30% increase in CPU over the 
basic Eulerian solution. 
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a ) 


b) 


c ) 


Figure 7.4: Vorticity contours at channel mid-height, a) initial, b) final with Eulerian 
scheme and c) final with Eulerian/Lagrangian scheme (increment = 2.0). 
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») b) c) 

Figure 7.5: Contours in S (entropy related function) at channel mid-height, a) initial, 
b) final with Eulerian scheme and c) final with Eulerian/Lagrangian scheme (increment 
= 0.003). 
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Figure 7.6: Velocity vectors at channel mid-height, a) initial, b) final with Eulerian 
scheme and c) final with Eulerian/Lagrangian scheme. 
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Figure 7.7: Velocity profiles across vortex for exact, Eulerian and Eulerian/Lagrangian 
solutions. 



Figure 7.8: Pressure profiles across vortex for exact, Eulerian and Eulerian/Lagrangian 
solutions. 


115 






Figure 7.9: Pressure coefficient as a function of the convection distance along the chan- 
nel, a) Eulerian solution, b) Eulerian/Lagrangian solution and c) exact solution. 
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Chapter 8 


Preservation of a Turbulent Inlet Velocity 

Profile in a Pipe 


The preservation of a turbulent inlet velocity profile in a straight pipe is analogous 
to the preservation of the Lamb vortex in a straight channel described in Chapter 7. 
Both the Lamb vortex and the turbulent velocity profile are characteristics of a viscous 
flow. The exact inviscid solution of the flow problem states that the velocity profile 
is not a function of the convection distance. As in the case of the Lamb vortex, the 
interest here is to see how these viscous flow features represented though the velocity 
profile are affected by the numerical diffusion of the Eulerian scheme. The turbulent 
velocity profile is provided as the inlet boundary condition of the calculation. As the 
viscous effects are not taken into account in the Euler equations, each flow particle 
should conserve its speed as it convects downstream. When using an Eulerian solver, 
the changes in the velocity profile are the result of numerical diffusion. The Lagrangian 
correction technique is then applied in order to recover the exact velocity profile at 
any station along the pipe. This case bears similarity with the convection of a fully 
developed turbulent flow in that the velocity profile is not a function of the convection 
distance. 

This test case is a preamble to the computation of secondary flows in bent pipes 
Chapter 11 where the source for the secondary flow creation is introduced as a non- 
uniform inlet stagnation pressure (obtained through a non-uniform inlet velocity pro- 
file), and where the Euler equations are used to predict the subsequent secondary flow 


117 




Figure 8.1: Straight circular pipe computational grid (161x25 nodes), 
development. 

Since the flow is incompressible and steady, the Euler equations are modified by 
the artificial compressibility method as described in Section 2.2. The Eulenan and the 
Eulerian/Lagrangian solutions are computed using the same grid, shown in Figure 8.1, 
with 161 nodes on an axial cross-section and 25 cross-sections equally spaced along the 
pipe length. Since the flow is symmetric, the computations are performed on only half 
of the pipe. The ratio of pipe length to pipe radius is 5.52. The inlet cross-section shows 
the unstructured grid with an O-type grid near the pipe wall and an H-type grid near 
the center of the pipe. This type of combined grid is chosen for the pipe geometries over 
the standard O-grid which presents a singularity at the pipe center. 


In the case of swirling flow in pipes or secondary flow calculations, the Lagrangian 
solution te chni que, using an upstream integration of streamlines and markers placed at 
the center of each cell, is preferred over the downstream integration technique, since the 
vorticity is diffused in the whole domain. 


The streamlines are also traced backward in the present test case, with two integra- 
tion steps per cell and are recomputed each 25 iterations of the Eulerian solver. The 
correction procedure is applied only once every Eulerian iteration and the corrections are 
multiplied by a factor 1 /4 in order to limit the perturbations imposed on the Eulerian 
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solution. 


A fully developed turbulent velocity profile is assumed at the pipe inlet cross-section. 
The inlet velocity profile corresponds to a Reynolds number of 

Re = — = 2.48 X 10 5 (8.1) 

V 

based on a mass-flow averaged velocity u of 2.31 m/s (chosen here as the reference veloc- 
ity) and on the pipe diameter d of 0.1023 m. The kinematic viscosity v is 0.952 x 10 -6 
m 2 /s. The pipe radius R is chosen as the reference length. The inlet velocity profile is 
determined by the ‘universal velocity distribution law’ for smooth pipes and very large 
Reynolds numbers given in [62] 


« + = = 5.75 log 10 (j/ + ) + 5.5, 

u 

where u * is the friction velocity defined as 



and is a Reynolds number based on u* as 

y + = — , 

v 

where y is the distance from the wall of the pipe. 


( 8 . 2 ) 


(8.3) 


(8.4) 


Since the flow computation is inviscid, the velocity at the wall must take a finite value 
called the ‘slip velocity’. Imposing a zero velocity value at the wall would tend to produce 
non-physical reverse flow under any small positive pressure gradient perturbation. In 
order to define the slip velocity, Prandtl’s universal law of friction for smooth pipes 
given in [62] is used to provide a relationship between the coefficient of friction A 

A = 8 (£) ! (8.5) 

and the Reynolds number as 

4= = 2.0 log 10 {ReV\) - 0.8, (8.6) 
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Figure 8.2: Universal velocity distribution law for smooth pipes. 


For a Reynolds number of 2.48 X 10 5 , this relation gives A = 0.01503. The shear stress 
at the wall r w and the friction velocity u* are now found by 


t w 



(8.7) 


and Equation (8.3). The slip velocity is chosen as the velocity just outside of the viscous 
sub-layer at a y+ value of 30. 


Using the values for r/+ and u*, the slip velocity is found from the ‘universal velocity 
distribution law’ given in Equation (8.2). The resulting slip velocity u s u p is equal 
to 60.7% of the flow averaged velocity at a distance from the wall of 0.0286% of the 
pipe radius. The universal velocity distribution law for smooth pipes is represented in 
Figure 8.2 as well as the slip velocity location. 


In Figure 8.3, the velocity profiles on the exit cross-section are plotted as a function 
of the normalized radius r/R for the exact solution (inlet velocity profile), the Eule- 
rian solution and the Eulerian/Lagrangian solution. Due to the numerical diffusion 
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of the scheme, the sharp gradients near the pipe wall are smoothed out in the Eule- 
rian solution. With the Lagrangian correction scheme, this effect is canceled and the 
Eulerian/Lagrangian solution and the exact solution are within drawing accuracy. 

As for the unsteady flow case, the Lagrangian correction could be applied during 
the pseudo- unsteady convergence of the Eulerian solution. However, for the steady flow 
cases, it is more effective in terms of CPU reduction to apply the Lagrangian correction 
once the Eulerian solution has converged (or nearly converged) since the convergence of 
the Eulerian solution takes more iterations than the Lagrangian correction convergence. 

For this particular test case, the combined Eulerian/Lagrangian solution requires 
only ~ 60 iterations to converge (from the converged Eulerian solution), whereas the 
Eulerian solution requires ~ 1000 iterations to converge (from a constant flow initial 
solution) to a maximum residual of 1 x 10~ 5 . The Eulerian solution, based on the 
Lax-Wendroff scheme and the artificial compressibility concept, takes 274 x 10~ 6 
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seconds/iteration/grid node. In comparison, the combined Eulerian/Lagrangian scheme 
takes ~ 8.0 seconds/iteration for this particular frequency of trajectory integration and 
marker number. The CPU increase due to the Lagrangian correction for this test case 
is 'V 30% of the basic CPU required for the Eulerian solution. 
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Chapter 9 


Swirling Flow in a Pipe 


Another preamble case to the secondary flow computation in bent pipes of Chap- 
ter 11 is the analysis of a swirling flow in a straight circular pipe. A swirling flow in a 
straight pipe superimposed on a uniform axial velocity is chosen as a representation of 
the convection of the secondary flow in the straight section of pipe downstream of the 
bend exit cross-section. 

The model for the swirling flow and the Eulerian and Eulerian/Lagrangian solutions 
are presented in Section 9.1. In Section 9.2, the strength of the swirling flow is increased 
and the vorticity field is shown to concentrate by the phenomenon of vorticity gradient 
augmentation to a point where the vorticity field becomes inaccurately represented 
when using the fixed Eulerian or Lagrangian spatial discretization. As a result, the 
convergence of the combined Eulerian/Lagrangian solution is affected. The sources for 
the vorticity concentration are identified in Section 9.2.1 and a solution to the problem is 
proposed in Section 9.2.2 as the introduction of a pseudo-diffusion term in the Helmholtz 
equation. 


9.1 Swirling flow model and solution 


As in the case of the secondary flow, the swirling flow of this example is composed of 
two counter-rotating vortices and, since the solution is symmetrical, it can be performed 
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Figure 9.1: Straight circular pipe computational grids: coarse grid with 125 x 25 nodes 
and fine grid with 384 X 49 nodes. 


on a half-pipe geometry. The ratio of pipe length to pipe radius is L/R = 5.52. 


The calculations are performed on a coarse grid with 25 cross-sections equally spaced 
along the pipe length and 125 nodes on an axial cross-section and on a fine grid with 49 X 
384 nodes. Both grids are shown in Figure 9.1. The present model for the inlet velocity 
profile is characterized by high cross-flow velocities near the pipe walls and the half-pipe 
symmetry surface according to the results found in bent pipes in Chapter 11. The inlet 
cross-flow velocities are found from the Poisson equation relating the streamfunction 
to the axial vorticity 

V 2 * = -a> a , (9.1) 


where the streamfunction is defined by 


* = -Cy{ 1 - 


y 2 + ■z 2 x 
R 2 ' 


(9.2) 


R is the pipe radius and C is a constant characterizing the strength of the swirling flow. 
The cross-flow velocities v and w are found by 


v = 


dV 2Cyz 




w = — - 


= 4 -^)- 


dz R 2 ’ dy ~ V* R 2 

The resulting axial vorticity distribution is linear with respect to y 


(9.3) 


U?x = 


-8 Cy 
R 2 * 


(9.4) 


The first chosen cross-flow strength is such that the maximum inlet cross-flow velocity 
is 40% of the uniform convection velocity, i.e. C = 0.2. The inlet cross-flow velocities 
and axial vorticity distribution of the coarse grid solution are shown in Figure 9.2. Since 
the vorticity is obtained as a cell-averaged quantity, the minimum and maximum inlet 
axial vorticity levels are different on the coarse and the fine grids. Nevertheless, the 
differences in minimum and ma ximum axial vorticity levels between the two grids are 
less than 3% of the maximum axial vorticity. 


The Eulerian and the Eulerian/Lagrangian calculations are performed on both grids. 
In the combined scheme, the markers are initially placed at the center of each cell and 
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Figure 9.2: Case C = 0.2: Cross-flow velocity vectors and axial vorticity on inlet cross- 
section for coarse grid calculations. 

the streamlines are integrated upstream with approximately two steps per cell. The 
streamlines are recomputed every 50 iterations of the Eulerian solver. The correction 
step is applied only once at each Eulerian iteration and the corrections are multiplied by 
a factor 1/4 in order to limit the perturbations to the solution during the pseudo-time 
integration process. 

Because of the induced velocities of one vortex on the other and the presence of 
the pipe wall, the center of each vortex is moving in an helicoidal pattern resulting in 
an alteration of the initially linear vorticity distribution. Figure 9.3 presents the axial 
vorticity contours on three cross-sections along the pipe for the Eulerian solution and the 
Eulerian /Lagr an gi an solution on the coarse grid and the fine grid. As the distance along 
the pipe increases, the constant axial vorticity contours undergo a rotation (stations 
xjL — 0.08 and xjL = 0.52). At the symmetry surface, the axial vorticity must be 
zero since dwjdy — 0 and v = 0 in the symmetry surface. This condition forces the 
axial vorticity contours to concentrate near the symmetry surface (stations xjL = 0.52 
and xjL — 0.98). The effect of numerical diffusion near the sy mm etry surface spreads 
the gradient of vorticity on approximately two cells in the Eulerian solutions. With 
the Lagrangian solution, the gradient of vorticity is concentrated on approximately 
only one cell and does not appear in the plots since the solution for the vorticity is 
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drawn only until the last cell center near the symmetry. The effect of the numerical 
diffusion is clearly seen in the coarse grid Eulerian solution as a spurious curving of the 
axial vorticity contours near the pipe wall (inducing smaller velocities on the pipe wall 
and therefore smaller circulation around the pipe). Using the Lagrangian correction, 
the axial vorticity contours are in a very good agreement with the Eulerian solution 
on the fine grid, especially near the pipe wall. The fine grid Eulerian solution still 
exhibits a small numerical diffusion effect near the pipe wall. When using the LagTangian 
correction on the fine grid, this effect is eliminated. 


The circulation T around a closed curve is an integral value used to quantify the 
numerical diffusion as a function of the distance along the pipe. In a barotropic flow 
field where viscous effects are non-existent or can be neglected and if the forces acting 
on the fluid are conservative, Kelvin’s theorem states that the circulation around an 
arbitrary closed curve moving with the fluid should remain constant. Using the material 
derivative, this is expressed as 



( 9 . 5 ) 


The convection of a material curve is obtained by setting ‘convective’ markers in a closed 
curve pattern at some location in the flow field once the solution has reached its steady- 
state. The material curve is deformed as each marker convects downstream with the 
local flow. Since these markers are used only for particle tracing purposes, they do not 
induce any correction in the Eulerian flow field. The tracing of the ‘convective markers 
is similar to the technique used for the ‘corrective’ markers trajectories integration. 
However, pseudo-time steps can not be used since the time-steps used for the integration 
of the stre amlin es have to be identical for every convective marker. The circulation is 
found by integration around the material curve as 


r = j> v • dr \ (9*6) 

where the integral is obtained by trapezoidal integration over the total number of con- 
vective markers. The material curve can be traced downstream or upstream and addi- 
tional convective markers are added to the material curve if the distance between two 
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Figure 9.4: Initial and final circulation contours. 


consecutive markers becomes too large due to the divergence of the streamlines. 

The chosen closed curve is initially placed around the inlet cross-section and the 
convective markers are traced downstream. Figure 9.4 shows the position of the initial 
and final locations of the material curve. The curve is initially located near the pipe wall 
where the largest numerical diffusion occurs. Figure 9.5 shows the circulation around 
the closed curve as a function of the convection distance along the pipe (the average 
location over all the convective markers). The Eulerian solution on the coarse grid 
exhibits the largest circulation change. By using the combined Eulerian/Lagrangian 
scheme on the coarse grid, the changes in circulation are reduced by a factor ~ 4 near 
the pipe exit. The circulations for the fine grid Eulerian solution and the fine grid 
Eulerian/Lagrangian solution are comparable since the correction is much smaller than 
in the case of the coarse grid solution. In the case of the fine grid solutions, the amount 
of diffusion is too small to judge the effectiveness of the Lagrangian correction (the 
comparison would have to be performed on a longer pipe in order to get a more diffused 
Eulerian solution). 

The coarse grid Euler solution requires ~ 2000 iterations to converge (from a con- 
stant flow initial condition) to a maximum residual of ~ 1 x 10~ 5 . The Eulerian/Lagrangian 
solution takes only ~ 275 iterations to converge (starting from the converged Eulerian 
solution). The computations are performed on a Stardent GS-2000 in vector mode. The 
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Figure 9.5: Case C = 0.2: Circulation on a closed curve as a function of the distance 
along the pipe, a) coarse grid Eulerian solution, b) coarse grid Eulerian/Lagrangian 
solution, c)fine grid Eulerian solution, d) fine grid Eulerian/Lagrangian solution. 
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combined scheme takes 2.8 seconds/iteration/node on average (for the marker number 
and frequency of streamline integration described above). Thus, in the coarse grid case, 
the increase in CPU due to the introduction of the Lagrangian correction is ~ 33% 
of the basic Eulerian solution. In comparison, the fine grid Eulerian solution requires 
~ 16 times more CPU than the basic coarse grid Eulerian calculation (starting from a 
uniform flow solution). 


9.2 Vorticity gradient augmentation 


As the strength of the swirling flow increases, the gradients of vorticity shown in Fig- 
ure 9.3 near the exit of the pipe become stronger. The evolution of the scalar field in a 
flow with high vorticity concentration has been previously studied by Knio in [38]. 

The phenomenon of vorticity gradient augmentation is intensified when the La- 
grangian correction is applied. Indeed, due to numerical diffusion, the Eulerian solution 
tends to smooth out strong gradients, which is not the case when using the Lagrangian 
correction. 

The intensification of the vorticity gradient is illustrated in Figure 9.6 where the 
strength of the swirling flow has been increased by a factor 1.5 over the previous test 
case (i.e. C — 0.3). The axial vorticity contours are plotted at stations along the pipe 
in Figures 9.6a) and 9.6b) for the coarse grid Eulerian and Eulerian/Lagrangian solu- 
tions and in Figures 9.6c) and 9.6d) for the fine grid Eulerian and Eulerian/Lagrangian 
solutions, respectively. When using the Lagrangian correction on the coarse grid, the 
gradient of vorticity at the exit of the pipe is increased compared to the Euleri an so- 
lution alone. A strong vorticity gradient region of a few cells width is created in the 
Eulerian /Lagrangian solution separating a high vorticity region from a low vorticity re- 
gion. The same behavior is seen in the fine grid solutions, but the gradient of vorticity 
is more intense and spread over one cell only in the Eulerian/Lagrangian solution. 
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In Section 9.2.2, the strength of the swirling flow is further increased (to C — 0.6) 
and the width of the vorticity gradient region decreases. The vorticity field becomes 
poorly represented (either on the fixed size grid in the Eulerian solution or by markers 
located at cell centers in the Lagrangian scheme). In turn, the poor sampling rate of 
the vorticity field leads to inaccurate corrections of the Eulerian solution and affects the 
convergence of the combined Eulerian/ Lagrangian scheme. 

In Section 9.2.1, the strain field and the vorticity field are first identified as the 
sources for the vorticity gradient augmentation. In Section 9.2.2, a solution to the 
problem is proposed as the introduction of a pseudo- diffusion term in the Helmholtz 
equation. 


9.2.1 Sources for vorticity concentration 

In Figures 9.7a) and 9.7b), the convergence of three stre amlin es from inlet to exit is 
shown for the coarse grid Eulerian and Eulerian/ Lagrangian solutions in the case of the 
swirling flow strength corresponding to C = 0.3. The streamlines Eire chosen such as to 
end up in the strong vorticity gradient region at the exit of the pipe. Since for this test 
case the source-terms for the vorticity are small (as will be shown later), the vorticity 
carried by the material particles along the streamlines is essentially constant and equal 
to the inlet vorticity at the inlet location of the markers. Because of the merging of the 
streamlines, the inlet vorticity assigned to each of the three markers differs by a large 
amount. When the streamlines merge along the pipe, a strong gradient in vorticity is 
created. When using the Eulerian scheme, the vorticity along the stre amlin es does not 
remain constant because of numerical diffusion and the creation of the strong vorticity 
gradient is inhibited. When the Lagrangian correction is used, the diffusion is reduced, 
and consequently the solution presents stronger gradients of vorticity. This effect is 
accentuated when using the Lagrangian correction since, the three streamlines show a 
stronger convergence in the Eulerian/ Lagrangian solution than in the Eulerian solution 
alone due to the stronger swirl effect. 
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Figure 9.7: Case C — 0.3: Frontal view of 3 streamlines drawn from inlet to exit of the 
pipe, a) Eulerian solution, b) Eulerian/Lagrangian solution (i=inlet, e=exit). 


The creation of strong vorticity gradients along the pipe can be explained from 
the action of the strain and vorticity field by using the Helmholtz equation giving the 
convective change in vorticity for an incompressible flow as 

^ = (d5-V)tT. (9.7) 


This equation states that the evolution of the vorticity attached to a particle convecting 
in an incompressible flow is governed by the generation of tilting and stretching source- 
terms along the particle trajectory. If the right-hand side terms (source-terms for the 
vorticity) are small, each component of <2 will approximately behave as a non-diffusive 
quantity as 



(9.8) 


By taking the gradient of Equation (9.8), relations are found which describe the behavior 
of the vorticity gradient for a convecting particle in a steady flow. 


■D(Vu^) 

Dt 


(v • V)Vu> x = V 



~o 


du du> x 
dr dx 


dv du} x 
dr dy 


dw du> x 

~d¥~d7' 


(9.9) 


134 



(9.10) 


v)y u , - ^ Duj y du du) y dv dtj} y dw du y 
Dt y ^ Dt dr dx dr dy dr dz ’ 


— <$. v)Vw - v dwdiv z 

Dt ^ Dt dr dx dr dy dr dz 

~o 

Or 


(9.11) 


J(V^x) 

Dt 

gga) 

Dt 


J>(V^z) 

Dt 


-(Vw* • V)iT - Vuj x x w, 

(9.12) 

-(Voij, • V)u - Vca y x CS, 

(9.13) 

-(Vu> x • V)u - Vu) z x u;. 

(9.14) 


Hence, if a scalar quantity (here any component of vorticity) remains approximately 
constant along a streamline, its gradient is governed along a streamline by both the 
local strain field V v find the local vorticity (rotation) field u;. 


In the case of a swirling flow through a straight pipe, the main component of vorticity 
is in the axial direction (a:). The axial vorticity and the gradient of axial vorticity along 
a streamline can be written as a function of convection time t. Prom Equation (9.7) 

w *(0 ~ w x(0) = ((u> • V)u)dt = S Ur , (9.15) 

and from Equation 9.12 

Vu> x (t) - Vu> x (0) = J (-(Vu> x • V)tf- Vu, x u >)dt = Sv u , m . (9.16) 

These equations represent the increase in axial vorticity and arial vorticity gradient 
along a streamline. Next, the right-hand side terms of Equation (9.15) and of the z- 
component of Equation (9.16) are integrated along two particular streamlines from inlet 
to exit. For these two streamlines, Equation (9.16) is verified by comparing the increase 
in the z-derivative of the vorticity at the exit of the pipe to the vorticity gradient 
computed from the solution at the pipe exit. Moreover, the source-term for the axial 
vorticity is shown to be indeed small compared to the source-term for the axial vorticity 
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2 - derivative as declared by hypothesis. 


In Figure 9.8a) two streamlines are shown in front and side views. The axial vor- 
ticity derivative du x /dz computed from the Eulerian/Lagrangian solution is plotted on 
the pipe exit cross-section with indication of the location of the end-points of the two 
streamlines. Figure 9.8c) represents the values of the right-hand sides of Equation (9.15) 
(5 Wr ) and of the third component of Equation (9.16) (S dux / dz ) (i.e. the source-terms 
for the axial vorticity and the source-terms for the z-derivative of the axial vorticity) 
integrated on the two streamlines from the pipe inlet to the pipe exit. The source-terms 
for the axial vorticity are indeed small for both stre am l in es compared to axial vortic- 
ity magnitude, therefore justifying the use of Equation (9.8). The magnitude of the 
source- terms for the z-derivative of the axial vorticity is larger for both streamlines. 

These terms, taken at the exit of the pipe, correspond to the values of vorticity deriva- 
tive displayed on the pipe exit cross-section in Figure 9.8b) at the indicated end-point 
location of the two streamlines (the inlet values for du > x /dz is 0). Thus, Equation (9.16) 
is verified for the particular case of these two streamlines and show that the strain and 
rotation fields combination integrated on the right-hand side of the equation is indeed 
the mechanism by which the large vorticity gradients are created along the pipe. Fur- 
thermore, these vorticity gradients appear in the absence of strong source-terms for the 
vorticity. For this test case, the vorticity is essentially convected passively with the 
material particles. 


9.2.2 Introduction of a pseudo-diffusion term 

Examining Figures 9.7a) and 9.7b), it can be seen that the information between two 
markers at the inlet is lost because of the lack of grid resolution. In order to transport 
information from approximately one marker per cell at the inlet, the grid resolution in 
the cross-flow plane would have to be increased by a factor ~ 4 near the pipe exit, see 
Figure 9.7a). As the strength of the swirling flow increases, so does the required grid 
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Figure 9.8: Case C — 0.3: a) Front and side views of pipe with 2 stre amlin es, b) 
axial vorticity derivative du) x jdz on exit cross-section with end-points of streamlines 
(inc.=0.5), c) source-term for axial vorticity: and source-term for du> x jdz: S du , x / dz 

along the 2 streamlines. 
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resolution. However, because the representation of vorticity gradient augmentation is 
inviscid, as the grid resolution increases so do the gradients of vorticity as shown in 
Figure 9.6. 

When the strength of the swirling flow is increased by a factor 2, i.e. C ~ 0.6, the 
vorticity gradients are further intensified to a point where the high vorticity gradient 
region covers less than one cell. (This region separates a high vorticity from a low vor- 
ticity). This phenomenon is shown in Figure 9.9 where contours in axial vorticity are 
drawn for diverse stations along the pipe. The vorticity gradient is intensifying along 
the pipe. When using the Eulerian scheme alone, the numerical diffusion spreads the 
high gradient of vorticity region over a few cells, whereas when using the Lagrangian 
correction, wiggles in the high vorticity gradient region appear because the represen- 
tation of this region by markers placed at cell centers (or by state-vectors at the grid 
nodes) becomes inadequate near the pipe exit (as seen when superimposing the grid 
on the exit-cross-section in Figure 9.9b)). As mentioned before, the insufficient spatial 
representation of the flow destabilizes the Eulerian solution through the correction pro- 
cedure. For example, the inaccurate velocity corrections lead to erroneous perturbations 
of the axial velocity through the solution of the continuity equation. 

As shown before, the use of larger size grids does not help to overcome the problem 
since the gradient region concentrates. Therefore, the gradients have to be ‘controlled’ so 
as to be supported on the fixed size grid. The goal of the Lagrangian correction technique 
is still the reduction of the numerical diffusion of the Eulerian solver without resorting 
to larger size grids. However, in the present flow cases, because of the vorticity gradient 
augmentation phenomenon, the additional control of vorticity gradients is required in 
order to get a stable Eulerian/ Lagrangian solution. 

The present section shows that the introduction of a pseudo-diffusion term in the 
He lmh oltz equation allows to control the strength of the vorticity gradients in the flow. 
If the diffusion level introduced by the Helmholtz pseudo-diffusion is maintained lower 
than the numerical diffusion of the Eulerian scheme, the resulting solution will still show 
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a substantial improvement over the Eulerian solution alone. 


The pseudo- diffusion term is chosen such as to take the same form as the real diffusion 
term arising in a viscous fluid, hence the Helmholtz equation is now written as 


Du) , x . __ o 

— = (u . V)v + i/fV 2 u>, 


(9.17) 


where the coefficient is the Lagrangian pseudo-diffusion coefficient proportional to 
the mesh size h. The ratio of the right-hand side terms is of the order of the pseudo- 
Reynolds number Rei and is a measure of the tendency of the vortex lines to be frozen 
to the fluid compared to the tendency of the vorticity to diffuse in the flow. Using V 
and L to represent the reference velocity and length, then 


(w • V)u V L 

viV 2 u? vi 


~ Re.L . 


(9.18) 


The pseudo- diffusion term is added to the tilting/stretching term and integrated along 
the marker trajectory using a predictor-corrector scheme. The evaluation of the La- 
grangian pseudo- diffusion term in the Helmholtz equation is similar to the evaluation of 
the smoothing term in the Lax-Wendroff algorithm of Section 3.2. However, the vector 
Q is known at the cell centers instead of the grid nodes. Therefore, the average values 
are first defined at the cell nodes by averaging the values of l3 at the eight surrounding 
cell centers. The Lagrangian pseudo-diffusion term is then evaluated at a cell center 
by su mmin g the differences between the eight nodes average values and the cell center 
value. The vorticity Laplacian is then interpolated from the cell centers to the marker 
location for each Lagrangian time-step. 


Definition of vorticity at a wall node 

When defining the node average values of the wall nodes receive contribution from 
only four cell centers. The vorticity value at the wall nodes is determined through the 
following relations defining the vorticity gradients at a node by a surface integral over 
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Figure 9.10: Estimation of the vorticity at a wall node. 


a pseudo-cell based on the cell centers as sketched in Figure 9.10 


Vu> x V = 

/ Vu x dV = 
Jv 

Js u; x nd5, 

(9.19) 

Vu, v F = 

[ VUydV = 

/ U JyftdSy 

(9.20) 


Jv 

Js 

Vw 2 V = 

[ Vto z dV - 
Jv 

J u ) z ndS, 

(9.21) 


where V is the volume of the pseudo-cell. Using the nomenclature defined in Section 3.1, 
the vorticity gradient at node 1 of Figure 3.1 is defined as the sum from the contributions 
from the eight surrounding cells A to H. The contribution of cell A is given by 


Vw xl = 

1 

4Vj\ 


u x fidS , 

(9.22) 

Vu> yi = 

1 

4Vi . 

Ai,7 2 ,As 

u? y nd5, 

(9.23) 

Vu > zl — 

1 

4 Vi , 


v z ndS, 

(9.24) 


and so on for the contributions from cells B to F. The value of C3 W at the wall node is 
then found from the nearest node in the direction perpendicular to the wall as sketched 
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(9.25) 


in Figme 9.10 by using the values of uj and Vw for this node as 


... d(2 

U! w — U> + 


where 


d d dx d dy d dz 

dr\ dx dr] dy dr] + dz dr] 


Figure 9.9c) shows the axial vorticity contours for the same cross-sections as Fig- 
ure 9.9a) and 9.9b) for the Eulerian/Lagrangian solution with a Lagrangian pseudo- 
diffusion term added to the Helmholtz equation. The axial vorticity contours are 
smoother showing that the the addition of the pseudo-diffusion terms allows to con- 
trol the vorticity gradients. The resulting vorticity field is more accurately represented 
on the fixed-size grid and the Eulerian/Lagrangian scheme is stable. As designed, the 
Lagrangian pseudo-diffusion term is shown to be very effective at smoothing out the 
strong vorticity gradients in the flow without affecting the smooth vorticity regions. 
The numerical value of vi used for this test case is 0.05. 

The amount of diffusion introduced in the solution is quantified by looking at the 
changes in circulation around a closed curve moving with the flow. The circulation 
around the closed curve initially placed around the inlet cross-section is shown in Fig- 
ure 9.11 as a function of the convection distance along the pipe (measured as a mean 
value of all the convective markers positions) for the Eulerian solution and the Eule- 
rian/Lagrangian solution with the addition of the Lagrangian pseudo- diffusion term in 
the Helmholtz equation. In Figure 9.11, the change in circulation when using an Eule- 
rian/Lagrangian scheme and a Lagrangian pseudo-diffusion term is shown to be small 
compared to the change observed in the circulation using the basic Eulerian solution. 

The convective change in circulation around a closed curve due to the introduction 
of a diffusion term in the Helmholtz equation can also be expressed as a contour integral 
of the curl of the vorticity as given in Appendix F. 

As the strength of the swirling flow increases, the phenomenon of vorticity gradient 
augmentation is responsible for the inaccurate representation of the vorticity in the flow 
which is the cause for a destabilization of the combined Eulerian /Lagrangian scheme. 
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The introduction of a pseudo-smoothing term in the Helmholtz equation has been shown 
to successfully smooth out the strong vorticity gradients and allows one to obtain a stable 
Eulerian/Lagrangian solution. By comparing circulations around a closed converting 
material curve, the Eulerian/Lagrangian solution gives still a substantial improvement 
over the Eulerian solution alone. The computation of the pseudo-smoothing term results 
in a ~ 5% increase of CPU over the Eulerian/Lagrangian basic solution. 
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Chapter 10 


Constant Stagnation Pressure Flow in a 90 ° 

Bend 


The incompressible flow in a 90° bend of circular cross-section with constant stagna- 
tion pressure inlet conditions is examined. The constant stagnation pressure is obtained 
at the inlet through constant pressure and constant velocity conditions. The inlet ve- 
locity is set as part of the inlet boundary conditions. The inlet surface is located far 
enough from the bend so that the upstream influence of the bend is negligible and the 
constant pressure condition is also obtained. 

The exact solution for the flow field is at constant stagnation pressure since the 
stagnation pressure remains constant along a streamline and each streamline is charac- 
terized by the same stagnation pressure. In turn, the vorticity in the domain is zero 
everywhere because the flow is irrotational upstream of the inlet cross-section. From 
[4], if at any instant in time before t = 0 (time at which particles cross the inlet section) 
the flow is irrotational, then the flow remains irrotational for any subsequent instant. 
That is if all derivatives of Q are zero for any time previous to t = 0, then 

Duj D 2 Q D n u; 

~Dt = "DP = *" = ~dF = °’ ( 101 ) 

and if all derivatives are defined, Taylor’s theorem shows that the quantity tU vanishes 
for all subsequent instants in time. Alternatively, Hawthorne showed in [26] that the 
growth of streamwise vorticity is a function of the stagnation pressure field only, so that 
no streamwise vorticity can be created in a constant stagnation pressure field. 
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Because of truncation errors, the Eulerian solution does not exactly satisfy the con- 
stant stagnation pressure and zero vorticity conditions. By computing Eulerian solutions 
on a coarse and a fine grid, this test case provides for a calibration of local and global 
errors in stagnation pressure and vorticity. Also, by using the Eulerian/Lagrangian 
scheme, the errors in vorticity and stagnation pressure are shown to be less than for the 
Eulerian solution alone. This test case also provides a check on the numerical integra- 
tion of the vorticity source-terms (uJ- V)u. These must be driven towards zero (starting 
from an Eulerian solution with errors in vorticity and stagnation pressure) smce the 
exact solution for the flow is irrotational. The improvement due to the addition of 
the Lagrangian correction is compared to the improvement obtained when using the 
Eulerian scheme on a finer grid. 

Again, since the flow is symmetrical, the computations are performed only one half 
of the pipe. The geometry used for this test case is taken from the Enayet et al. [21] flow 
data and will also be used in Chapter 11 for secondary flow calculations. The geometry 
and grids are identical to the ones used for the accuracy study of Section 3.9. 

The pipe diameter is 0.048 m and the ratio of radius of curvature to pipe diameter 
is 2.8. The computational domain extends two diameters upstream of the bend inlet 
cross-section and two diameters downstream of the bend exit cross-section. As men- 
tioned before, the distance of two diameters upstream is required to satisfy the constant 
stagnation pressure condition at the pipe inlet cross-section. The calculations are per- 
formed on a coarse grid with 43 stations along the bend total length and 189 nodes 
in a cross-section. The fine grid is formed by 85 cross-sections with 713 nodes in a 
cross-section. The fine grid is composed of eight times as many cells as the coarse grid. 
The two grids are shown in Figure 10.1 in side and frontal views. The cross-sections at 
which the results are presented are indicated. Eulerian solutions are computed on both 
grids, whereas an Eulerian/Lagrangian solution is computed on the coarse grid only. 

The Lagrangian correction uses the upstream streamline integration scheme where a 
marker is placed at the center of every cell and the streamlines are integrated backwards 
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Figure 10.1: Coarse and fine grids front and side views (189 X 43 nodes and 713 X 85 
nodes) with particular cross-sections. 
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with approximately two steps per cell until they reach the inlet. The streamlines are 
recomputed every 50 iterations of the Eulerian solver. The correction step is applied 
once at each Eulerian iteration. An under-relaxation on the vorticity correction is used 
according to Equation (6.48) in order to limit the perturbations to the Eulerian solver. 
The under-relaxation factor Rf is 0.5%. 

Since the Lagrangian correction technique is based on vorticity corrections, the local 
errors are first shown in terms of vorticity errors. Contours in the z-component of vor- 
ticity are shown in Figure 10.3 for a cross-section located at 45° in the bend. The coarse 
grid Eulerian solution and the coarse Eulerian/Lagrangian solution are shown in Fig- 
ures 10.3a) and 10.3b). The fine grid Eulerian solution is represented in Figure 10.3c). 
Similarly, in Figures 10.4, 10.5 and 10.6, the contours of the z-component of vorticity 
are drawn for three cross-sections located at 90° (bend exit cross-section) and at 1 and 
2 diameters downstream of the bend exit cross-section (the z-component of vorticity is 
streamwise for these stations). The maximum amount of error in the cross-section is 
indicated in each case. The vorticity is normalized with respect to the group U{ n /R 
representing the inlet velocity divided by the pipe radius. 

The errors in vorticity are concentrated near the wall of the pipe. By compar- 
ing the three solutions, it is clear that the combined Eulerian/ Lagrangian scheme is 
very effective at reducing the errors in vorticity in the flow field. Indeed, in the Eule- 
rian/LagTangian solution, the error in the z-component of vorticity is reduced to nearly 
zero, whereas using the Eulerian scheme on a finer mesh leads to very little reduction 
of the error in vorticity. In the coarse grid Eulerian solution, the errors in vorticity are 
more diffused in the domain, but the level of error remains approximately the same for 
the coarse and the fine grid. 

The results are also calibrated in terms of the local error in stagnation pressure 
coefficient A C Po defined at any node by Equation (3.75). 

Figure 10.7 shows contours in A C Po on the 45° cross-section for the coarse grid Eule- 
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rian and Eulerian/Lagrangian solutions and the fine grid Eulerian solution. The 90°, Id 
and 2d cross-section results are represented in Figures 10.8, 10.9 and 10.10, respectively. 
The maximum error is indicated for each cross-section. The error in stagnation pressure 
is the largest at the inside of the bend near the 90° station where the pressure recovery 
occurs. By using the Eulerian/Lagrangian scheme, the error in stagnation pressure is 
reduced locally when compared to the reference case of the coarse grid Eulerian solution. 
The maximum error in stagnation pressure for the Eulerian/Lagrangian solution is even 
lower than for the fine grid Eulerian solution. Since the stagnation pressure is not trans- 
ported along the streamlines in the Lagrangian scheme, the correction of the stagnation 
pressure occurs indirectly through the vorticity correction (Crocco’s equation relates 
vorticity and stagnation pressure gradient) as mentioned in Section 4.2. The corrected 
velocities alter the flux balance around each cell in the Lax- Wendroff algorithm so that 
the pressure is corrected too. 

A global indicator of stagnation pressure losses is the L 2 norm of the stagnation 
pressure losses for the domain defined as e po by Equation (3.74). 

Even if the maximum local error in stagnation pressure is smaller with the Eule- 
rian/Lagrangian scheme than with the fine grid Eulerian solution, the errors are less 
diffused with the latter so that an average measure over the nodes like the L 2 norm 
gives a more unfavorable result for the Eulerian/Lagrangian solution. The L 2 norm 
of the stagnation pressure errors decreases by a factor ~ 2.3 when using the Eule- 
rian/Lagrangian solution, whereas the Eulerian solution on the finer grid leads to a 
reduction factor of ~ 4.0 (which is in accordance with the Lax- Wendroff scheme second- 
order accuracy). 

In an inviscid flow, by expressing the gradient of pressure in terms of the local 
coordinates system ( s,n,b ) where s, n and b stand for the streamwise, normal and 
binormal directions as shown in Figure 10.2, the gradient of pressure along the binormal 
direction can be shown to be zero [39]. This also implies that the norm of the velocity 
presents no dependence on the binormal direction since the flow is at constant stagnation 
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Figure 10.2: Local coordinate system (s,n,b). 

pressure for this test case. Also, because the secondary vorticity is zero (as is the case for 
the other components of vorticity), the streamlines do not present any torsion, and cross 
the grid cross-sections perpendicularly. The contours in velocity norm are represented 
in Figures 10.11 and 10.12 for the 45° and 90° stations, respectively. These plots 
provide an indication of the vorticity correction near the pipe wall. The exact solution 
presents no dependence on the binormal direction (here the direction perpendicular to 
the symmetry surface, i.e. the y-direction). The reduction of the spurious deformation 
of the contours of velocity norm is larger when using the Eulerian/Lagrangian scheme 
them when using the Eulerian scheme on the fine grid. 

The Eulerian/Lagrangian scheme is clearly more efficient than the fine grid Eulerian 
solution at reducing the errors in vorticity in the flow. Indeed, while the vorticity 
errors are more concentrated when using a finer grid, the level of vorticity remains 
approximately the same in the coarse or the fine grid Eulerian solution. Also, it has been 
shown that the vorticity correction of the Lagrangian scheme results in a more accurate 
solution of the velocity field than the fine grid Eulerian solution. The m aximum error 
in stagnation pressure in the field is lower when using the Eulerian/Lagrangian scheme 
than the fine grid Eulerian scheme. However, since the stagnation pressure errors are 
more diffused in the flow with the Eulerian/Lagrangian scheme, the reduction in the X 2 
norm of the stagnation pressure losses is lower than with the Eulerian solution on the 
fine grid (a factor 2.3 compared to a factor 4.0). 
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The Eulerian/Lagrangian solution is substantially less expensive than the fine grid 
Eulerian solution. The addition of the Lagrangian correction leads to an increase of 
~ 75% of the basic coarse Eulerian solution. In comparison the fine grid Eulerian solu- 
tion requires 16 times more CPU than the coarse grid Eulerian solution (the Eulerian 
solution on the coarse grid requires ~ 2000 iterations to converge to a maximum resid- 
ual of ~ 1 x 10 5 (the initial flow is interpolated from a coarser grid solution), the 
Eulerian/Lagrangian solution requires ~ 390 iterations to converge starting from the 
converged Eulerian solution and takes ~ 10.6 seconds/iteration in average). 
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Figure 10.3: Contours of z-component of vorticity at 45° station, a) Eulenan solution, 
b) Eulerian/Lagrangian solution, c) Eulerian solution on fine grid (increment = 0.05). 



Figure 10.4: Contours of streamwise vorticity at 90° station, a) Eulerian solution, b) 
Eulerian/Lagrangian solution, c) Eulerian solution on fine grid (increment = 0.05). 
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a) 

(AC pl ) mM = 2.5% 


b) 

(A C P0 ) max = 0.9% 


c ) 

(AC Po ) max = 1.0% 


Figure 10.7: Contours for the local error in stagnation pressure coefficient A C Po at 45° 
station, a) Eulerian solution, b) Eulerian/Lagrangian solution, c) Eulerian solution on 
fine grid (increment = 0.005). 




(A C PQ ) max = 10.6% (A C P0 ) max = 2.7% 



c ) 

(AC Po ) mox = 4.9% 


Figure 10.8: Contours for the local error in stagnation pressure coefficient A(7 po at 90° 
station, a) Eulerian solution, b) Eulerian/Lagrangian solution, c) Eulerian solution on 
fine grid (increment = 0.005). 


153 



a ) b) c ) 

{& C Po)rnax = 4.6% (AC Po ) max = 0.8% (A C Po ) max = 1.3% 

Figure 10.9: Contours for the local error in stagnation pressure coefficient AC Po at Id 
station, a) Eulerian solution, b) Eulerian/Lagrangian solution, c) Eulerian solution on 
fine grid (increment = 0.005). 



a ) b) C ) 

(A^po)mai = 4.4% (AC Po ) max = 0.9% (AC' Po ) max = 1.1% 

Figure 10.10: Contours for the local error in stagnation pressure coefficient AC Po at 2d 
station, a) Eulerian solution, b) Eulerian/Lagrangian solution, c) Eulerian solution on 
fine grid (increment = 0.005). 
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a) b) c ) 

Figure 10.11: Contours in velocity norm at 45° station, a) Eulerian solution, b) Eule- 
rian/Lagrangian solution, c) Eulerian solution on fme grid (increment = 0.01). 


0.92 



a) b) c) 

Figure 10.12: Contours in velocity norm at 90° station, a) Eulerian solution, b) Eule- 
rian/ Lagrangian solution, c) Eulerian solution on fine grid (increment = 0.01). 
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Chapter 11 


Secondary Flow in Bent Pipes 


The secondary flow phenomenon results from the existence of a vorticity component 
in the ‘throughflow’ direction. As the flow moves along a bend, each fluid particle is 
subjected to a centrifugal force acting in a direction of a line going through the center of 
curvature of the streamline. If the velocity profile is non-uniform, the centrifugal force 
acting on faster moving particles will be larger than the one acting on slower moving 
particles. This induces a motion of the particles in a cross-flow plane with non-uniform 
cross-flow velocities. For example, a fully- developed pipe flow profile presenting higher 
velocities at the pipe center and lower velocities near the pipe wall will lead to secondary 
flow when taken through a bent pipe. Under the action of the centrifuged forces, the 
slow moving fluid is pushed inwards whereas the fast moving fluid moves towards the 
bend outer wall. Because of continuity, the fast moving fluid displacement towards the 
outer wall forces slow moving fluid to convect along the pipe wall towards the pipe inner 
wall. A circular motion constituting the secondary flow phenomenon is then established 
in a cross-flow plane. This is illustrated in Figure 11.1 taken from [19]. The effect of the 
secondary flow in bent pipes of circular cross-section is then to displace the high velocity 
regions towards the outer wall. But if the inlet boundary- layer is thin, the main core flow 
will behave approximately as a free-vortex flow through the bend and the displacement 
of particles in cross-flow planes will be small compared to the displacement along the 
streamwise direction. 

As a result of the Helmholtz equation, the secondary flow in bent pipes can be 
explained in terms of the turning, stretching and diffusion of the vorticity attached to a 
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Figure 11.1: Secondary flow in a bent pipe due to non-uniform streamwise velocity. 

given fluid particle. For the flow in a bent pipe, the vorticity may be initially present in 
the inlet velocity profile. In the case where the flow is fully developed at the inlet of a 
pipe of circular cross-section, the inlet vortex lines form concentric rings. If the diffusion 
term is ignored, the secondary vorticity is due only to the turning and stretching of the 
legs of the inlet vortex lines which are parallel to the pipe symmetry-plane, as two 
particles on the inner and outer radii convect through the bend with different speeds. 
This is illustrated in Figure 11.2 where a few vortex lines are drawn from an Eulerian 
solution of a flow in a 90° bend. At the inlet of the computational domain, the flow is 
fully developed and the vortex lines present no throughflow component. Near the bend 
inlet, the few vortex lines drawn near the wall surface indicate the presence of secondary 
vorticity due to tilting and stretching. The phenomenon intensifies as the flow moves 
down the bend as shown by the vortex lines drawn near the bend exit. 


The creation of the secondary flow in the bend has been described so far as an 
‘inviscid’ flow phenomenon. The additional effect of the developing boundary-layers or of 
the wall log region in the case of a fully developed flow is to counteract the development 
of the inviscid secondary flow as explained in [21]. The boundary-layer on the inner wall 
of the pipe is subjected to a favorable pressure gradient, so that it grows slowly. On the 
outer wall the adverse pressure gradient causes the boundary-layer to thicken. The result 
of taking the boundary-layers into account is then to displace the high velocity regions 
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more towards the inner wall, an effect inverse to the prediction of the inviscid secondary 
flow theory. If the bend exit is under constant pressure conditions (for example due to 
a straight section or an exit at atmospheric pressure) the streamwise pressure gradients 
have to be reversed near the exit in order to readjust to the downstream pressure. The 
evolution of the boundary-layer is therefore complex to predict. 

The secondary flow phenomenon is a feature common to many flow problems and 
has been the object of extensive studies. Both approximate solutions and numerical 
computations have been tested against experiments. Berger et al. [10] review article 
provides an extensive list of references for the flow in curved pipes. Hawthorne [28] 
studied the applicability of secondary flow analyses for the solution of internal flow 
cases. Ackeret [2] also described peculiarities linked to the internal flow behavior. 

A well-known approximate solution for the secondary vorticity in a bend or a cascade 
of airfoils is the Squire and Winter relation [67], which links the secondary vorticity 
generation to the bend angle and the inlet vorticity. Rowe [59] compared experiments 
and computations based on the Squire and Winter formulation for a 180° bend and 
found reasonable agreement up to an angle of about 75°. Detra [19] also derived an 
approximate solution procedure based on s m all perturbations of the throughflow in 
good agreement with experiments in pipes of 21° and 42° bend angle. However, the 
pipe radius has to be small compared to the radius of curvature for this approximate 
solution to be valid. 

More recently, Lakshminarayana [39] derived generalized expressions for the sec- 
ondary vorticity using intrinsic coordinates. This work has been extended by Hawthorne 
[29] for stratified fluids in rotating systems. 

An analytical formulation based on a streamlike function formulation of the inviscid 
flow in a curved duct was proposed by Abdallah [1] and tested on a 90° duct of rectan- 
gular cross-section. The turning of the velocity contours is predicted in good agreement 
with experiments. Briley [12] investigated three-dimensional viscous flows with large 
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secondary vorticity by deriving a system of approximating equations for flows where a 
primary direction is present, but in which the transverse velocity components can be 
large. Also, an extensive series of laminar and turbulent Navier- Stokes solutions in bent 
pipes of circular cross-section and square ducts have been compared with experiments. 
A few of these are [35], [34], [33], [74], [79], [32], [52]. 


11.1 Motivation 


The present work uses the Euler equations for the prediction of secondary flows in bent 
pipes. The motivation for the use of an inviscid formulation as reported by Hawthorne 
in [28] is based on the fact that large regions of the flow, although non-uniform, may 
be assumed frictionless since the influence of the walls extends only gradually inwards 
as the flow passes around the bend. Also, the presence of viscosity is not required 
for the generation of secondary vorticity even if the initial presence of vorticity at the 
inlet of the pipe is due to a viscous effect. The fundamental behavior of the flow can 
be predicted using the Euler equations since the development of the secondary flow is 
mainly the result of an inviscid process of vortex stretching and tilting. This applies for 
high Reynolds number flows since the ratio of time scales between the transverse viscous 
momentum diffusion and the convection time scale is proportional to the inverse of the 
Reynolds number for laminar flows. In turbulent flows, the ratio is roughly constant 
but still small. 

Also, the Euler equations are not restricted by the small shear assumption used in 
the approximate solutions of [67] and [19] where the distortion of the vorticity by the 
secondary vorticity is neglected as a second-order effect. This is particularly important 
here, where high turning angles and high Reynolds numbers are of interest. 

An underlying motivation is the study of the numerical diffusion effects on the so- 
lution. By comparing Eulerian solutions on coarse and fine grids, it is shown that the 
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numerical diffusion inherent to the Eulerian solver has an important effect on the so- 
lution both in terms of strength and location of the secondary flow. The fact that the 
general effect of the numerical diffusion is shown to be qualitatively similar to the real 
viscous effects is a misleading result since the solution remains very sensitive to the 
amount of numerical diffusion. Therefore, very fine grids are required for secondary 
flow computations in order to obtain a solution in which the numerical diffusion effect 
is small. By computing the circulation around a closed curve moving with the fluid, it 
is shown that even the fine grid Eulerian solution still suffers from numerical diffusion, 
indicating the limitation of standard Eulerian solvers for the computation of secondary 
flows. The motivation for the introduction of the Lagrangian correction technique is to 
find an alternative to expensive large size grids Eulerian calculations of secondary flow 
phenomena. With the addition of the Lagrangian correction technique, the amount of 
numerical diffusion observed in the Eulerian solution is shown to be reduced for a given 
grid size when compared to the standard Eulerian solution, or equivalently the grid 
requirements are reduced to converge to a solution where the effects of the numerical 
diffusion are small. 

The reduction of the spurious numerical diffusion to an unimportant level allows 
one to address the question of how much the secondary flow is influenced by the real 
fluid viscosity by comparing the truly inviscid numerical solution with the existing 
experimental data. 

Earlier, Navier- Stokes solutions have been reported to suffer from numerical diffu- 
sion more than from the turbulence modeling uncertainties [35]. This was due to the 
coarse discretisation imposed by large memory size requirements. Today, the turbulence 
modeling together with near wall boundary conditions is likely to be the predominant 
uncertainty factor in the viscous solution of the secondary flow in bends as seen in [32]. 
In comparison, Eulerian/Lagrangian computations offer a computationally cheaper ap- 
proach and, as seen in the result sections, still provide with a description of the basic 
flow behavior. However, the limitations of the present method are listed below. 
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11.2 Limitations 


As mentioned earlier, the viscosity effects tend to reduce the development of the sec- 
ondary flow. Therefore, inviscid numerical solutions tend to overpredict the secondary 
flow generation, a result consistent with the findings of this study. The effect of the 
growing boundary-layers is to narrow the channel and change the pressure distribution. 
The shear stress for a laminar boundary- layer subjected to a streamline curvature can 
be written as a function of 6/R c where 8 is the boundary-layer thickness and R c is the 
radius of curvature of the pipe. For a turbulent boundary- layer, however, the effect 
of the streamline curvature on the boundary-layer thickness are an order of magnitude 
larger [11]. In the turbulent flow case, the curvature of the bend acts as a flow destabi- 
lizer near the outer wall whereas the flow is stabilized near the pipe inner wall [2]. At 
the same time, the fluid is transported from inner to outer wall by the secondary flow 
motion. 

Another limitation for the use of the Euler equations is the requirement that a 
slip velocity has to be defined at the wall for the inlet velocity profile (if not reverse 
flow would immediately occur on the outside of the bend under the influence of the 
adverse pressure gradient). However, the problem is not singular to the Euler equations 
solutions, since the requirement of a ‘cut-off’ velocity at the wadi also arises in the 
approximate solutions methods because the small shear assumption becomes invalid 
near the wall. Squire and Winter [67] used a ‘cut-off’ velocity of 45% of the mean velocity 
and Detra [19] tried both 65% and 80% cut-off velocities on the inlet velocity profile. 
As reported in [28], the limitation is important since the secondary vorticity depends 
on the incoming velocity profile vorticity. Hawthorne [27] found a cut-off velocity by 
canceling the overpredicting effect of the secondary vorticity of Squire and Winter by 
a large negative vorticity at the wall. However, this analysis is limited to small bend 
angles. 

The present approach to the slip velocity definition for a turbulent inlet mean velocity 
profile is to neglect the laminar sublayer. The slip velocity is then chosen as the velocity 
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just outside the laminar sublayer as mentioned in Chapter 8 where this approach has 
been used for the preservation of a turbulent inlet velocity profile in a pipe. The model 
is valid for turbulent boundary-layers where the thickness of the boundary-layer is much 
larger than the laminar sublayer. The present choice of = 30 as the upper limit of 
the laminar sublayer is arbitrary. However, this value is consistently used to determine 
the slip velocity for the two flow cases computed in this Section. Still, the best way 
to determine this slip velocity would be to calibrate the t/+ value by fitting numerical 
computations to measurements. This was not attempted here, since a wider study would 
be required in order to determine the best value of y+ valid for a range of flow cases. 
Instead, as mentioned in Section 11.1, the present chapter deals with the reduction of 
the spurious numerical diffusion phenomenon. 

Denton [18] uses a similar concept in his inviscid calculations of turbomachinery in 
order to introduce viscous effects in the momentum equation by the use of a distributed 
body force. The wall shear stress is obtained by neglecting the displacement thickness 
of the laminar sublayer, that is the surface streamline lies at the edge of the laminar 
sublayer. A loss function is derived as a power law distribution from the wall into the 
field. A value of t/+ of 10 for the edge of the laminar sublayer is recommended by [18] 
only if there are enough nodes to describe accurately the boundary-layer. In practice, 
a value between 10 and 40 produced the best fit with measurements. 

A first attempt to counteract the overprediction of the secondary-flow is performed 
in this study by introducing a simple ‘wall correction 1 where the velocity is adjusted 
at every node on the wall, using the universal log-law distribution between the nodes 
adjacent to the wall and the nodes on the wall. This method is described in more detail 
in Section 11.4.3. 

The comparison of numerical solutions with experiments is, however, suffering from 
a lack of experimental data on the cross-flow itself for the cases treated here. Refraction 
of the laser beam at the water-plexiglass interface prevented the measurements of cross- 
section velocity components in the pipe of circular cross-section of the Enayet data 
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set [21]. The data available is essentially composed of streamwise velocity contours for 
diverse cross- sections along the bend supplemented by pressure measurements along the 
wall for diverse angles around the bend. The absence of experimental measurements of 
the cross-flow (necessary to quantify directly the strength and position of the secondary 
flow) makes the comparison with the calculations incomplete. However, calculations 
are helpful since they indicate a correlation between streamwise velocity contours and 
position and strength of the cross-flow regions. 


11.3 Outline 


The first pipe geometry is taken from the Enayet et al. [21] data set while the second 
corresponds to a bend tested at the Gas Turbine Laboratory of MIT. 


First, coarse and fine grid Euler solutions of the flow are compared to indicate 
the effect of numerical diffusion on the strength and position of the secondary flow. 
The Lagrangian correction technique is then introduced in order to reduce the level of 
numerical diffusion, allowing the identification of the real viscous effects by comparison 
with experiments. The numerical diffusion reduction is demonstrated by computing the 
circulation around bend cross-sections along the pipe and by tracing the circulation on 
a closed curve moving with the fluid. Then, results are shown using the wall correction 
method which attempts to counteract the overprediction of the secondary flow. Pressure 

measurements along the bend are also compared with experiments for the Enayet et al. 
geometry. 
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11.4 Enayet 90° bend case 


The secondary flow development for an incompressible flow (water) through a 90° bend 
of constant circular cross-section is investigated. The results axe compared to experimen- 
tal data taken from the Enayet et al. [21] data set providing Laser-Doppler Velocimeter 
measurements of throughflow or ‘streamwise’ velocities at different cross-stations along 
the bend and pressure measurements along the wall for diverse angles around the bend. 

The geometry is identical to the one used for the constant stagnation pressure flow 
in Chapter 10. The diameter of the pipe is 0.048 m. and the ratio of radius of curvature 
to pipe diameter is 2.8. The computations were performed on a coarse and a fine grid 
with 320 X 36 nodes and 1223 X 71 nodes, respectively. The grids front and side views 
are represented on Figure 11.3. The grid cross-sections are spaced evenly in the stream- 
wise direction along the bend centerline. The test geometry extends three diameters 
downstream of the bend exit cross-section. In order to accelerate the calculations, the 
computational grid reaches only two diameters downstream. It has been verified by 
numerical experiments that the computational results are not affected by applying the 
exit boundary condition at such a reduced distance from the exit of the bend. The grid 
extends 0.58 diameters upstream of the bend inlet cross-section where the inlet velocity 
is known from measurements. Eulerian solutions are performed on both grids, whereas 
the Eulerian/Lagrangian solution is computed on the coarse grid only. 

A particularity of this secondary flow problem is that the vorticity (either inlet 
boundary-layer vorticity or secondary vorticity) presents concentrated regions but is 
also present in a more diffused form in a large portion of the flowfield. This implies 
that the Lagrangian correction technique has to operate on a large extent of the flow 
therefore requiring a large amount of markers. The option to place a marker at the 
center of each cell has been chosen here for simplicity. The upstream integration of 
the streamlines is found to be more successful than the downstream integration since it 
implies an even distribution of the markers in the flow. The correction step is applied 
only once at each Eulerian iteration. Both an under-relaxation of the correction (R f = 


165 




d side views (320 x 36 nodes and 1223 X 71 


0.1%) and a multiplication of the corrections by a factor 1/4 are used in order to limit 
the perturbations to the Eulerian solution. The streamlines are recomputed each 50 
iterations of the Eulerian solver. In order to prevent the formation of strong gradients 
in vorticity, a Lagrangian pseudo-diffusion coefficient of 2% is used. 

The Enayet data set provides streamwise velocity measurements for a turbulent flow 
case at a Reynolds number of 43000 and a flow averaged velocity u of 0.92 m/s. The 
fluid kinematic viscosity is 0.804 X 10 — ® m^/s. The stations at which the experimental 
data are available are located at 30°, 60° and 75° along the bend and at 1 diameter 
downstream of the bend exit section as indicated in Figure 11.3. 


11.4.1 Inlet velocity profile definition 

The measured inlet velocities on horizontal and vertical traverses at a distance of 0.58 
diameter upstream of the bend inlet section are splined in the circumferential and radial 
directions to give an inlet boundary condition for the calculation. The slight asymmetry 
in the inlet velocity measurements has been averaged out for this calculation. The inlet 
velocity profile is not a fully developed profile due to the proximity of the bend but the 
assumption of a fully developed turbulent flow is used here for the determination of the 
slip velocity at the wall. The procedure is identical to the one described in Chapter 8. 
The coefficient of friction, shear stress and friction velocity take the values 

A = 0.022, t w = 0.00233, Ur = 0.0482. (11.1) 

Using a value of y + = 30 or u + = 14 at the edge of the laminar sublayer, the slip 
velocity u s u p is 

u slip = u + u T = 0.675, or ^ = 0.730. (11.2) 

The distance from the wall at which the last inlet velocity measurement is taken is 
equal to 10% of the pipe radius. Since the first grid point out from from the wall is 
located in between the last experimental value and the wall, its velocity has to be fitted 
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Figure 11.4: Inlet velocity profile at grid nodes. 


on the velocity profile. The velocity profile is represented at the grid nodes location in 
Figure 11.4. 

In order to be able to compare the Eulerian solution on a coarse and on a fine grid, 
both calculations must have the same inlet boundary condition. When discretizing 
an analytical velocity profile onto a coarse grid, the inlet vorticity is lower than the 
vorticity obtained on a fine grid. This problem is illustrated in Figure 11.5 for the case 
of an inlet velocity profile of boundary-layer type. To ensure that the vorticity values 
between grid nodes is identical between the two cases, the velocity profile is modified 
when using the fine grid according to the arrows in Figure 11.5. The inlet vorticity 
value is particularly important here since it is used as a component of the state- vector 
in the Eulerian/Lagrangian solution and also as a measure of the numerical diffusion of 
the solution when comparing Eulerian or Eulerian/Lagrangian results on a coarse grid 
to an Eulerian result on a fine grid. 
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Figure 11.5: Velocity profile defined on coarse and fine grids and adjustment of the 
velocity profile for the fine grid case. 
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11.4.2 Enayet case 90° bend: Eulerian and Eulerian/Lagrangian re- 
sults 

This section presents the secondary flow solution for the Eulerian solution computed on 
the coarse grid, the Eulerian/ Lagrangian solution on the coarse grid and an Eulerian 
solution on the fine grid. The three solutions are compared with experiments. 

The streamwise velocity contours and the cross-flow velocity vectors are shown for 
the four stations along the bend on Figures 11.6, 11.7, 11.8 for the Eulerian solution 
on the coarse grid, the Eulerian/Lagrangian solution and the Eulerian solution on the 
fine grid, respectively. The experimental streamwise velocity contours taken from [21] 
are represented on Figure 11.9. 

The cross-flow velocity vectors indicate a weak secondary flow region is present 
along the pipe wall at the 30° station. On further stations, this region intensifies, moves 
towards the inside of the bend, and then moves along the pipe symmetry surface as 
the vortical region is entrained due to the the proximity of the wall and the symmetry 
surface. 

Generally, the Eulerian solution on the fine grid and the Eulerian/Lagrangian solu- 
tion on the coarse grid lead to a stronger secondary flow than the Eulerian solution on 
the coarse grid. The presence of the pipe wall and of the counter-rotating vortex entrains 
the vorticity around the pipe wall and then along the symmetry surface. Because of the 
difference in the secondary flow intensity, the entrainment speed and the location of the 
vortical regions is also different in the three cases. Therefore, the numerical diffusion 
affects both the strength and the location of the secondary flow regions. Using the Eule- 
rian solution on the coarse or the fine grid results in very different vortical flow location 
at the Id station as can be seen on Figure 11.6 and 11.8. With the Eulerian solution 
on the coarse grid, the secondary flow region is weaker due to numerical diffusion and 
the vortical region moves more slowly along the wall and the symmetry plane. When 
using the Eulerian/Lagrangian solution on the coarse grid, the vortical regions move 
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faster than with the Eulerian solution alone and axe closer to the Eulerian solution on 
the fine grid predictions showing that the combined scheme is successfully at reducing 
diffusion errors. 

Not only is the strength and location of the vorticity field affected by the numerical 
diffusion for the Eulerian solution on the coarse grid, but the structure of the vortical 
field is also affected. With the Eulerian solution on the fine grid, the main secondary 
vortex present on the Id station is accompanied by two smaller counter-rotating vortices 
located near the inside of the bend. With the Eulerian solution on the coarse grid, the 
two small vortices are smeared out due to numerical diffusion. With the Lagrangian 
correction, the two small vortices are recognizable at the Id station indicating that the 
structure of the secondary flow is closer to the one predicted by the fine grid Eulerian 
solution. 

For each of the three solutions, the low streamwise velocity regions are associated 
with the strong cross-flow regions so that the location of the cross-flow in the experiments 
can be estimated by looking at the experimental velocity contours on Figure 11.9. For 
both the Eulerian/ Lagrangian solution on the coarse grid and the Eulerian solution 
on the fine grid, the predicted speed of entrainment of the cross-flow seems to be too 
large indicating that the neglected viscous effects near the pipe wall are important. The 
Eulerian/Lagrangian solution shows a good agreement with the Eulerian solution on 
the fine grid. In particular, Figure 11.7 shows how the gradients in axial and cross-flow 
velocities near the pipe wadi and the symmetry plane are better represented when using 
the Eulerian/Lagrangian solution. 


Nevertheless, in the Eulerian solution on the coarse grid, the location of the cross- 
flow is in better agreement with the experiments because the numerical diffusion plays 
a role in decreasing the strength of the cross-flow and hence its entrainment speed. The 
numerical diffusion helps the solution when comparing the solution with experiments in 
terms of cross-flow strength and location. However, this is a misleading result since the 
numerical diffusion and the real viscous diffusion are two distinct phenomena. Indeed the 
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Eulerian solution shows strongly diffused streamwise velocity contours, a result which 
does not agree with the experimented streamwise velocity. Therefore the numerical 
diffusion of the flow is a spurious computational effect which has to be minimized. 

Figures 11.10, 11.11 and 11.12 show the streamwise vorticity contours on the four 
stations for the Eulerian solution on the coarse grid, the Eulerian/Lagrangian solution 
and the Eulerian solution on the fine grid, respectively. The effect of the numerical 
smoothing is clearly visible in Figure 11.10 (Eulerian solution on the coarse grid) by 
the curving of the streamwise vorticity contours near the wall of the pipe. This effect 
is reduced for the solution on the fine grid in Figure 11.12. The use of the combined 
Eulerian/Lagrangian scheme also results in the correction of the diffusion effect near 
the pipe wall as seen in Figure 11.11. 

Nevertheless, the Eulerian/Lagrangian solution on the coarse grid still suffers from 
lack of grid resolution when compared to the Eulerian solution on the fine grid. The 
effectiveness of the Eulerian/Lagrangian technique is limited because the flow features 
are too small to be accurately captured on the coarse grid. If the Eulerian solution on 
the fine mesh is taken as the representation of the inviscid non-diffused solution one 
can see that the secondary vorticity gradients to be captured are very high and that 
the effectiveness of a solution on a grid twice as coarse is limited by a grid resolution 
issue. This is the case when the formed secondary vorticity zone moves along the pipe 
symmetry line. 

Figure 11.13 shows the streamlines emerging from the inlet boundary-layer near 
the wall and wrapping around the vortex for the Eulerian solution on the coarse grid 
(a), the Eulerian/Lagrangian solution (b) and the Eulerian solution on the fine grid 
(c). In the Eulerian solution on the coarse grid, because of the numerical diffusion, the 
vortex is formed further downstream compared to the two other solutions and the rate 
of rotation around the vortex is lower. The Eulerian /Lagrangian solution predicts the 
vortex formation even earlier than the Eulerian solution on the fine grid. This is because 
in the first part of the pipe, the Lagrangian correction scheme is more effective than the 
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fine grid Eulerian solution in preventing numerical diffusion. This is not true, however, 
in the second part of the pipe where the vorticity is very concentrated and where the 
Lagrangian correction effectiveness is limited by the grid resolution issue mentioned 
earlier. 

Figure 11.14 compares the pressure contours for the three solutions on the half- 
pipe symmetry surface. The negative streamwise pressure gradient on the inner wall is 
followed by a pressure recovery. The trace of the secondary vortex on the symmetry 
surface is clearly visible in the fine grid Eulerian solution. Up to 2/3 of the pipe length, 
the Eulerian /Lagrangian solution shows a good agreement with the fine grid Eulerian 

solution. 

Figure 11.15 shows the circulation around pipe cross-sections as a function of their 
axial distance downstream s (measured on the pipe axis and normalized by the pipe 
radius R) for the Eulerian solution on the coarse grid (a), the Eulerian/Lagrangian 
solution (b), and the Eulerian solution on the fine grid (c). The condition of constant 
circulation around pipe cross-sections does not apply since the contours defined around 
pipe cross-sections do not correspond to the convection of a material curve initially 
placed around the inlet cross-section. However, this allows for an integral value com- 
parison between coarse and fine grid calculations. The circulation (normalized by the 
mass-flow averaged velocity u and the pipe radius) increases from the inlet until approx- 
imately the 70° station due to the creation of the secondary flow. With the solution 
on the coarse grid, the circulation increases more slowly due to the strong presence of 
numerical diffusion damping the strength of the secondary flow. The circulation com- 
puted with the Eulerian/Lagrangian scheme, however, agrees better with the fine grid 
Eulerian solution. The departure between the Eulerian/Lagrangian solution and the 
Eulerian solution on the fine grid increases after the 60° station, approximately. This is 
an indication of the numerical diffusion occurring on the finer mesh. The upstream in- 
fluence of the bend is shown by the increase in circulation for the straight pipe segment 
placed before the bend. 
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An interesting circulation contour to assess the numerical diffusion of the secondary 
flow is a material curve surrounding the exit cross-section and convecting ‘upstream’ 
until approximately the 90° cross-section location. The value of the circulation should 
remain constant while the material curve is convecting. Again, the circulations for the 
three types of calculations are plotted in Figure 11.16. The symbol s now stands for 
the average convection distance from the exit cross-section. The numerical diffusion 
is seen as an increase of the circulation when the material curve convects ‘upstream 5 . 
Also because of the numerical diffusion the Eulerian solutions start from lower values of 
circulation and the two circulations increase as the material curves convect upstream. 
In comparison, the Eulerian/Lagrangian solution shows little change in the circulation. 
Still the curve for the Eulerian solution on the fine grid shows a substantial amount 
of loss in circulation. This integral value plot indicates that very fine grid Eulerian 
calculations are indeed required for these particular vortical flows in order to get a truly 
non-diffused solution. 

The convergence towards a solution where the numerical diffusion is small results 
in an increase of the entrainement speed of the secondary vortex. The predicted speed 
of the vortical regions for the Eulerian/Lagrangian solution and the Eulerian solution 
on the fine grid are too large when comparing the computed results to the streamwise 
velocity measurements. The difference is accounted for by the effects of the real fluid 
viscosity effects. An attempt to take these effects into account is described in the next 
section. 
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Figure 11.6: Stream wise velocity contours and cross-flow velocity vectors for four sta 
tions along the bend using Eulerian scheme on coarse grid. 
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Figure 11.8: Streamwise velocity contours and cross-flow velocity vectors for four sta- 
tions along the bend using Eulerian scheme on fine grid. 














Figure 11.13: Streamlines emerging from the near wall region at the inlet of the pipe 
and forming the secondary flow: a) Eulerian solution, b) Eulerian/Lagrangian solution, 
c) Eulerian solution on fine grid. 


v/ 



Figure 11.14: Pressure contours on half-pipe symmetry surface (increment = 0.02): a) 
Eulerian solution, b) Eulerian/Lagrangian solution, c) Eulerian solution on fine grid. 
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Figure 11 . 16 : Circulation around a closed convecting curve: a) Eulerian solution, b) 
Eulerian/Lagrangian solution, c) Eulerian solution on fine grid. 
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11.4.3 Law of the wall correction 


As mentioned in the previous section, the streamwise velocity measurements are in- 
dicative that the inviscid solution of the flow predicts a speed of entrainment of the 
secondary flow which is too large. While the Eulerian solution predicts the fundamental 
flow features, the effect of the fluid viscosity is to decrease substantially the strength of 
the secondary flow and therefore its entrainment speed. The effects of the fluid viscosity 
were briefly described in Chapter 11. This section presents an attempt at accounting 
the shear effects near the wall of the pipe as part of the flow solution. 


In order to account for the viscous phenomenon occurring at the wall a simple Taw 
of the wall’ is imposed on each wall node. The log-layer region of the boundary-layer 
is in equilibrium, i.e. it does not depend on the pressure gradient imposed upon it. 
Under the assumption that the first grid point out from the wall is embedded in the 
log-layer, the wall velocity depends only of the velocity value at the first neighboring 
node since the log-layer is described by a universal distribution law. The correction at 
a wall node consists then in scaling the velocity by using the same log-law relation than 
for the inlet velocity profile described in Chapter 8. For each wall node, the magnitude 
of the velocity at the first grid point out from the wall is used to determine the value 
of the friction velocity u* by implicitly solving the universal velocity distribution law 
given in Equation (8.2). Each component of the slip velocity at the wall is then scaled 
so that the resulting velocity magnitude is equal to the velocity magnitude at the end 
of the laminar sublayer at a u+ value of 14 (this value is identical to the one used for 
the inlet velocity profile definition). 

Figure 11.17 shows the streamwise velocity contours and the cross-flow velocity vec- 
tors on the four stations along the pipe for the Eulerian/Lagrangian solution on the 
coarse grid with the addition of the simple wall velocity correction. When using the 
Eulerian/Lagrangian scheme, the markers in the cells near the pipe walls are omitted 
since the wall velocity correction procedure already determines a vorticity value m the 

cell. 
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The streamwise velocity contours are now in better agreement with the experimental 
values drawn in Figure 11.9. The addition of the wall correction leads to a reduction 
of cross-flow velocities near the pipe wall and symmetry surface. This means that the 
mechanism of near wall diffusion is indeed a factor to take into account when simulating 
the flow field. The agreement is particularly good at the 75° station. After this station, 
the main vortex departs from the near wall region and travels along the symmetry 
surface. At the Id station, the vortex seems to be located higher than in the experiments 
showing the absence of numerical diffusion modeling in the interior of the flow field. 

In general, the agreement is not as good when the viscous effects are not taken 
mto account as shown in the previous section. However, the simple ‘law of the wall 
correction’ suffers from the absence of viscous diffusion out from the wall regions. Also, 
the assessment of the quality of the solutions proves to be difficult since the present 
experiments do not provide any direct information on the cross-flow itself. 
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Figure 11.17: Streamwise velocity contours and cross-flow velocity vectors for four sta 
tions along the bend using Eulerian/Lagrangian scheme on coarse grid and wall correc 
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The Enayet data set also provides pressure measurements along the wall of the pipe 
for 5 angles around the bend from the bend inside at 0° to the bend outside at 180° 
and for 6 stations along the bend located at 15°, 30°, 60°, 75° and at 1 and 2 diameters 
downstream of the bend exit cross-section. 


Figures 11.18 shows the computed pressure along the wall for angles around the 
bend at 0°,45°, 90°, 135° and 180° for the Eulerian solution on the coarse grid (a) and 
the Eulenan/Lagrangian solution (b). The pressure is non-dimensionalized by the use 
of a reference pressure p re f (computed as the average pressure on the inlet cross-section) 
and a dynamic head based on the flow averaged velocity u. The abscissa represents the 
curvilinear coordinate along the bend measured on the pipe axis. The measurements are 
indicated by different symbols for each angle around the bend. The measurements show 
a pressure loss of approximately 0.3 pu 2 from inlet to exit. Since the Eulerian solver does 
not take into account any physical loss mechanism, the pressure is entirely recovered 
at a distance of 2 diameters downstream of the bend exit. The agreement between 
measurements and calculations is only qualitative. Both calculations overpredict the 
pressure for the 0°,45° and 90° angles. The agreement is better at 135° but on the 
bend inside (at 180°) the pressure recovery is much larger than in the experiments. 
In the experiments, the pressure recovery region on the bend inside corresponds to an 
increase in the boundary-layer thickness neglected in the Eulerian calculations. 


Figure 11.19 shows the computed pressure along the wall when the wall velocity 
correction is used in the Eulerian/Lagrangian solution. The loss mechanism at the 
wall has generally little effect on the pressure distribution. While the agreement has 
improved from a solution without wall correction, the coincidence with experiments 
remains only qualitative. 
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11.5 GTL 90° bend case 


The second test geometry has been tested at the Gas Turbine Laboratory of MIT [6]. 
Streainwise velocity measurements are provided at a station located at 1.61 diameters 
downstream of the bend exit cross-section. 

The diameter of the pipe is 0.1023 m. and the ratio of radius of curvature to pipe 
diameter is 1.61. The computations were performed on a coarse and a fine grid with 
320 X 51 nodes and 1223 X 101 nodes, respectively. The grids front and side views are 
represented on Figure 11.20. The grid cross-sections are spaced evenly in the streamwise 
direction along the bend centerline. The test geometry extends 1.61 diameters upstream 
and downstream of the bend exit cross-section. The upstream distance is required for 
the influence of the bend to be negligible at the inlet cross-section since the inlet velocity 
profile is specified as a fully developed flow. 

The flow Reynolds number is 3.32 x 10 5 and the flow averaged velocity is 2.61 m/s. 
Using the procedure described in Chapter 8, the coefficient of friction, wall shear stress 
and friction velocity take the values 

A = 0.0142, t w = 0.0121, u T = 0.1101. (11-3) 

Using a value of y + = 30 or u + = 14 at the edge of the laminar sublayer, the slip 
velocity u,u p is 

u,u p = u + Ur = 1.541, or ^ = 0.589. (11A) 

The upstream integration of the streamlines is used in the Eulerian/Lagrangian 
technique. The correction step is applied only once at each Eulerian iteration. Both an 
under-relaxation of the correction {R f = 0.1%) and a multiplication of the corrections 
by a factor 1/4 are used in order to limit the perturbations to the Eulerian solution. 
The streamlines are recomputed every 50 iterations of the Eulerian solver. In order to 
prevent the formation of strong gradients in vorticity, a Lagrangian pseudo-diffusion 


191 




coefficient of 5% is used. 


11.5.1 GTL case 90° bend: Eulerian and Eulerian/Lagrangian results 

The streamwise velocity contours are shown for the station located at 1.61 diameter 
downstream of the bend exit in Figure 11.21 for the Eulerian solution on the coarse grid, 
the Eulerian/Lagrangian solution and the Eulerian solution on the fine grid, respectively. 
The experimental streamwise velocity contours are represented on Figure 11.21 d). 

The three numerical solutions again predict an entrainment speed of the secondary 
flow which is too large. This is especially true of the fine grid Eulerian solution and 
the Eulerian/Lagrangian solution. The latter is overpredicting the secondary flow en- 
trainment by the largest amount indicating that this solution is the one with the least 

numerical diffusion. 
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11.6 Conclusions for the secondary flow in bends 

The use of the Euler equations for the solution of the secondary flow in bent pipes has 
been motivated by the fact that viscous effects should appear in only limited regions 
of the flow. Indeed, the Euler solution of the problem has been shown to predict the 
fundamental features of the flow for the two bend cases treated here. However, the 
solution is largely dependent on the amount of numerical diffusion, even with the fine 
grid solution. By using the Lagrangian correction technique, the numerical diffusion on 
the coarse grid was reduced to a level below the one observed on the fine grid. Based on 
this ‘true inviscid’ solution, the assessment of the fluid viscosity effects can be assessed 
by comparing with experiments. One of the limitations of the method is the definition 
of a slip velocity at the wall. 

Consistently with the findings of [19], [67], and [28], the inviscid solution of the flow 
has been shown to overpredict the development of the secondary flow. A law of the wall 
correction has been implemented to take the near wall viscous effects into account. The 
resulting streamwise velocity contours and cross-flow position compared better with the 
measurements. However, the pressure distribution along the wall was only in qualitative 
agreement with measurements. The concordance with the experiments also deteriorated 
when the secondary flow moved away from the wall since no viscous effects are taken into 
account apart from the wall regions. The law of the wall correction is also dependent 
on the chosen value of y + at the edge of the laminar sublayer. Clearly, this area of 
research requires more study. For example, combining an expression like the power law 
of Denton to a calibration of the value of y + over several experiments would yield more 
promising results. 
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Chapter 12 


Weston Wing Case 


The external flow over a three-dimensional wing tested by Weston [78] at the Langley 
Research Center is the object of this section. The wing geometry is characterized by a 
rectangular untwisted planform, a NACA0012 cross-section, a semi-spam to chord ratio 
b/c of 3 and a body of revolution tip. 

The Euler equations have been previously used by many authors to describe the flow 
around wings [45, 55, 36, 50, 51, 30, 9, 77]. These equations provide realistic solutions 
for these flows, since vorticity is captured as part of the solution and since the dynamics 
of the wake roll-up and convection is essentially inviscid. It is generally accepted that 
the artificial dissipation is the cause for the flow separation at sharp tr ailin g edges. 
Moreover, as reported by Roberts [56], separation has been observed also on rounded 
wing tips, a phenomenon believed to be linked to grid resolution. 

The ffeestream Mach number and angle of attack are 0.1425 and 8°, respectively. 
The experimental data consists of wake and pressure coefficient measurements. The 
flow around the identical geometry has been computed by Roberts [56], on two grids of 
different topology, with a pressure coefficient on the wing in good agreement with the 
experiments. However, his calculations showed a spurious diffusion of the tip vortex 
behind the trailing edge when compared with experiments. This effect was identified as 
a numerical diffusion artifact. 

Because of the low Mach number, the approximation of incompressible flow is used 
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for the present computations. The Eulerian solution of the flow is first used for a 
validation of the incompressible flow solver by comparing the calculation to experiments. 
The Lagrangian correction is then applied here in order to counteract the numerical 
diffusion observed in the trailing vortex region. 

The grid with a C-H structure is generated as described in Appendix A. The grid 
has a C structure at each spanwise cross-section and the wing surface is discretized by a 
H grid structure as shown in Figure 12.1. Figure 12.2 shows a detail of the definition of 
the wing tip body of revolution surface. Each streamwise grid surface in the wake has a 
H structure and presents a strong clustering in the wing tip vortex region as shown in 
Figure 12.3. This type of grid in the wake has been selected instead of the 0-0 type of 
grid used by Roberts because of the higher resolution of the wake. The mesh is slanted 
in the wake so that the high clustering region follows approximately the trailing vortex 
upward movement behind the trailing edge. The grid extends 2.3 chords downstream 
of the tr ailing edge and 2.5 chords away from the wing tip in the spanwise direction. 
The minimum distance from the wing to the grid outer surface is 2.5 chords. Since the 
farfield boundary conditions are based on the normal velocity component through the 
boundary, the external surface of the grid was inclined at an angle with respect to the 
freestream in order to get a non-zero normal velocity component. 


The farfield boundary conditions use the theory developed in Section 3.4. On the 
farfield boundary, the freestream velocity at an incidence angle of 8° is corrected by the 
amount of velocity induced by the horseshoe vortex system. 


The computed pressure coefficient on the wing surface is presented at five span- 
wise locations and compared with the experimental values of [78] in Figure 12.4. The 
computed pressure coefficient is defined as 


r „_ P ~P°° 
yvi’ 


( 12 . 1 ) 


where the subscript refers to the freestream values, is interpolated linearly from the 
grid nodes to each chosen spanwise location. Computed results and experiments are in 
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Figure 12.1: Weston wing grid 101 x 26 x 17 nodes with C-H structure shown by 2 mesh 
surfaces. 
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Figure 12.2: Weston wing grid detail of leading edge-tip region. 



Figure 12.3: Exit surface with clustering near trailing vortex region and symmetry mesh 
surface showing wake surface angle behind trailing edge. 
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good agreement. In this calculation as in Roberts’, the predicted leading-edge suction 
peak is higher near the wing root than over the rest of the wing. Roberts related this 
effect to the considerable tunnel flow angularity near the wing root. Also, the present 
solution gives lower leading-edge suction peaks near the wing root and higher near 
the wing tip compared to Roberts’ solution. This can be attributed to a coarser grid 
resolution near the wing root and a better resolution of the tip geometry. Generally, 
Roberts’ calculation and this calculation both lead to an underprediction of the suction 
peak over most of the wing, when comparing with the experiments. Again, this can 
be related to the flow angularity in the tunnel. However, the underprediction near the 
wing root is more pronounced in the present calculation since the wing resolution is 
approximately four times coarser than in the solution of Roberts. The grid resolution 
on the wing was traded for more grid resolution in the wake when the C-H grid type 
was selected instead of the 0-0 type. 

A singular grid line extends from the wing tip to the outer surface corresponding to 
a locally lower solution accuracy and accounting partially for the larger discrepancies 
between experiments and calculations in this region. The deterioration of the computed 
solution near the wing tip was also reported by Roberts and the local solution in the 
tip region was shown to be very sensitive to the details of the grid and wing geometry. 

The Lagrangian correction technique is applied on the Eulerian solution by placing a 
marker at the center of the cells in the vicinity of tr ailin g vortex and by integrating the 
streamlines backward towards the tr ailin g edge. An under-relaxation factor Rf = 1% 
and a multiplication of the corrections by a factor 1% are used in order to limit the 
perturbations to the Eulerian solver. The correction step is applied once each Eulerian 
iteration. The integration of the streamlines is repeated each 5 iterations of the Eulerian 
solver. In contrast with the computation of the flow in pipes of Chapter 11 where the 
vorticity is spread over a large portion of the computational domain, the Lagrangian 
correction technique shows its potential in this case since the markers need to be located 
in only a small region of the flow and the CPU time allotted to the Lagrangian correction 
remains correspondingly low. 
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Figure 12.5: Initial location of markers in the trailing vortex region. 


Figure 12.5 shows the initial patch of markers before the integration of the stream- 
lines. Each marker is indicated by a dot. The markers are placed within a ‘conical 1 
region roughly accounting for the trailing vortex expansion downstream of the trailing 
edge. Because the i mm ediate region behind the trailing edge corresponds to very high 
gradients, attempts to correct this zone with the Lagrangian method (where vorticity 
is constant over a cell) were unsuccessful. The vortex is spread over very few cells and 
discretizing it by markers placed at cell centers is very inaccurate. Thus, the markers 
were placed only from 0.2 chords downstream of the trailing edge. 


The experimental survey of the wake provides axial vorticity, pressure coefficient 
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and stagnation pressure coefficient defined as 


Cpo — 


PO POoo 

IpOT’ 


( 12 . 2 ) 


for two stations in the wake located at 0.5 chords and 2.0 chords downstream of the 
trailing edge, respectively. Figure 12.6a) shows axial vorticity contours predicted by 
the Eulerian solver for the two stations downstream of the trailing edge. Figure 12.6b) 
shows the axial vorticity contours for the same stations when the Lagrangian correction 
is used. The experimental values of axial vorticity are reported in Figure 12.7. The 
maximum computed level of vorticity is generally lower than the experimental values, 
indicating that a finer grid around the wing and in the wake would be required to get 
the vorticity level encountered in the experiments. However, the reduction of the wake 
diffusion using the Lagrangian technique is still relevant. The experimental data show 
a decrease of vorticity by a factor 1.2 between the two stations, whereas the maximum 
vorticity decreases by a factor 1.5 between the two stations with the Eulerian solution 
alone, an effect of numerical diffusion. With the Lagrangian correction, the maximum 
vorticity level is higher at 0.5 chords downstream of the trailing edge indicating that 
the numerical diffusion is reduced. This level remains approximately unchanged in the 
wake. At the 2 chords station, the corrected and uncorrected maximum vorticity levels 
vary by a factor 2, indicating that the Lagrangian correction is capable of handling 
large corrections of vorticity. Also the vortex core is tighter when using the Lagrangian 
correction. 


Figure 12.8 shows the computed pressure coefficient contours on the two stations 
for the Eulerian solution and the Eulerian/Lagrangian solution. The experimental val- 
ues are shown in Figure 12.9. The diffusion of the vortex is clearly seen at the 2 
chords cross-section with the Eulerian solution. The experimental minimum pressure 
coefficient remains approximately unchanged between the two stations. With the La- 
grangian correction, the minimum pressure coefficient is closer to the experimental value 
and decreases between the two stations, an effect which will be explained later. The 
computed stagnation pressure coefficient is shown in Figure 12.10 and the experimental 
values are reported in Figure 12.11. Again, the minimum stagnation pressure coefficient 
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increases in the wake with the Eulerian solution due to numerical diffusion The same 
decreasing behavior seen for the pressure is also valid for the stagnation pressure in the 
Eulerian/Lagrangian solution. 


The maximum vorticity in the wake is traced as a function of the distance down- 
stream of the trailing edge for the Eulerian and the Eulerian/Lagrangian solutions in 
Figure 12.12. The diffusion of the Eulerian solution is seen between the trailing edge 
(z/c = 1.) and the exit cross-section ( x/c = 3). With the Lagrangian correction, the 
maximum vorticity remains approximately constant. The roughness of the curve is due 
to the vortex converting through the grid. As mentioned earlier, the vorticity correc- 
tion is not performed near the trailing edge but begins at 0.2 chords downstream of the 
trailing edge. 

In general the present solution shows pressure coefficient and stagnation pressure co- 
efficients in better agreement with experiments than Robert’s solution. This is believed 
to be due to the higher grid resolution in the vortex vicinity used in the present work. 
Roberts Eulerian solution also showed excessive numerical diffusion in the wake with 
an increase in stagnation pressure coefficient by a factor 3 between the two stations. 

The minim um pressure coefficient traced as a function of the distance downstream 
of the trailing edge is shown in Figure 12.13. The diffusion of the vortex leads to a 
increase of the minimu m value of the pressure coefficient with the Eulerian solution. 
The Lagrangian correction leads to a very different behavior. When the correction 
begins at 0.2 chords downstream of the trailing edge, the minimum pressure coefficient 
first decreases because more vorticity is entrained in the vortex and the larger vorticity 
region leads to a lower pressure at the vortex center. After 1 chord downstream of the 
trailing edge (z/c = 2.0), most of the trailing wake vorticity has been pulled into the 
vortex, and hence the m i nimu m pressure in the vortex re mains approximately constant. 


The numerical diffusion of a trailing tip vortex is shown by an Eulerian calculation. 
The Lagrangian correction is successful at preventing the vortex diffusion downstream of 
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the trailing edge. However, the loss of vorticity is the highest in the trailing edge region 
and the lack of resolution impeded the use of the Lagrangian correction at this loca- 
tion. A more effective correction procedure could be obtained by using the Lagrangian 
correction technique in conjunction with an adaptive grid refinement procedure, so that 
the vortex resolution remains approximately the same even near the trailing edge. The 
location of regions to be adapted would be indicated simply by the presence of the 
Lagrangian markers. Searching for ‘features’ in the Eulerian solution woud likely be 
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Figure 12.6: Axial vorticity in wake for stations located at 0.5 and 2.0 chords down- 
stream of the trailing edge, a) Eulerian solution, b) Eulerian/Lagrangian solution (inc. 
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Figure 12.7: Experimental axial vorticity in wake for stations located at 0.5 and 2.0 
chords downstream of the trailing edge. 
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Figure 12.10: Stagnation pressure coefficient in wake for stations located at 0.5 and 2.0 
chords downstream of the tra il i n g edge, a) Eulerian solution, h) Eulerian/Lagrangian 
solution (inc. = 0.05). 
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Figure 12.12: Maximum axial vorticity in wake as a function of the distance behind the 
trailing edge, a) Eulerian solution, b) Eulerian/Lagrangian solution. 



Figure 12.13: Minimum pressure coefficient in wake as a function of the distance behind 
the trailing edge, a) Eulerian solution, b) Eulerian/Lagrangian solution. 






Chapter 13 


Conclusions 


13.1 Summary 

A new approach is proposed as the coupling of an Eulerian and a Lagrangian solu- 
tion procedures for the reduction of numerical diffusion encountered in finite- difference 
time- marching Eulerian calculations. The motivation behind this work is the efficient 
numerical treatment of flow non-homogeneities, such as vortex wakes, embedded in an 
otherwise smooth background flow field. The coupling of the Eulerian and the La- 
grangian solution techniques is intended to enhance the poor vorticity and entropy 
capturing capabilities of standard Eulerian solvers. 

The first part of this thesis deals with the numerical methodology used for the so- 
lution of the Euler and Lagrange equations. The Euler equations are solved using a 
Lax-Wendroff algorithm on an unstructured grid, and are used for compressible and 
incompressible (through the artificial compressibility concept) flow situations. The 
numerical smoothing is based on a second- difference for compressible flows and on a 
fourth- difference (second-order accurate on distorted grids) formulation for incompress- 
ible flows. The far-field, wall and symmetry boundary conditions are described in both 
cases. A numerical study is performed to prove the second-order accuracy of the scheme. 

The Lagrangian solution method, based on particle markers, is then described. The 
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integration in time of the position, vorticity and entropy attached to the markers flowing 
at the local speed through the Eulerian grid is performed using a predictor/corrector 
scheme. The flow quantities required at a marker location are tri-linearly interpolated 
from the cell containing the marker to the marker location. The comparative advantages 
of two strategies of trajectory integration and/or marker positioning in the flow are 
described for representative flow situations. The coupling sequence between Eulerian 
and Lagrangian solutions is completed (and this is the key point of the procedure) by 
a ‘correction step’ in which the Lagrangian markers provide information on how to 
‘correct’ the spurious numerical diffusion of the Eulerian solution. As mentioned in this 
work, the correction step takes place only locally, each marker influencing only the cell 
where it is located. The correction procedure is described as a ‘vorticity correction 
where the vorticity attached to the marker is used to alter the velocity components at 
the nodes of a cell. Additionally, an ‘entropy correction’ is described for the compressible 
flow cases. The corrections are implemented as iterative procedures. The convergence 
speed of the vorticity correction process is reported. 

In all the test cases, the improvement due to the combined Eulerian/Lagrangian 
solver is tested by comparison to the standard Eulerian solution and to the Eulerian 
solution on a finer grid (where the numerical diffusion effects are lower). All combined 
Eulerian/Lagrangian solutions are computed on the basic coarse Eulerian grid. Also, in 
some flow cases, a comparison with experiments is performed. The CPU requirements 
for the combined scheme are also compared to the CPU needed for coarse and fine grid 
Eulerian calculations. 

The first test case is a compressible unsteady calculation, namely a Lamb vortex con- 
vection in a straight channel. The numerical diffusion effects are assessed by comparing 
the vortex solutions at the beg innin g and at the end of the channel. 

The preservation of a turbulent inlet velocity profile in a straight pipe shows the 
numerical diffusion occurring near the pipe wall (where the gradients are high) with the 
standard Eulerian solver. The use of the combined Eulerian/Lagrangian solver results, 
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instead, in a solution indistinguishable from the exact solution. 


The calculation of a swirling flow through a straight pipe, superimposed on a uniform 
axial velocity, is used as a preamble to a secondary flow calculation in a bent pipe. Again, 
numerical diffusion is identified in the Eulerian solution, especially near the pipe wall. 
By comparison with an Eulerian solution on a finer grid, the combined scheme is shown 
to be successful at reducing numerical diffusion errors. The circulation around a closed 
convecting curve (which should remain constant in an incompressible inviscid flow) 
is also used as an integral measure of numerical diffusion and shows the substantial 
improvement obtained with the combined scheme. A basic problem of inviscid flow 
calculations is also identified as a vorticity gradient augmentation, encountered in this 
swirling flow case, but also present when computing secondary flows in bent pipes. This 
concentration leads to a poor representation of the vorticity on the fixed size grid and 
destabilizes the combined scheme. The source for the phenomenon is identified as a 
‘vorticity convection process’ along streamlines, and appears in the absence of strong 
source-terms for the vorticity. The intensification process is shown to worsen when using 
the combined scheme because the reduction of numerical diffusion leads to a gradient 
definition over fewer cells. A solution is proposed as the introduction of a Lagrangian 
pseudo- diffusion term, similar in form to the true viscous diffusion term, introduced in 
the right-hand-side of the Helmholtz equation. This procedure is shown to result in 
a solution accurately described on the fixed grid. The change in circulation around a 
closed convecting curve is shown to be still much lower for the combined scheme with 
the Lagrangian pseudo- diffusion term than for the basic Eulerian solution. 

A constant stagnation pressure flow in a 90° bend is investigated. The errors in 
vorticity and stagnation pressure are reduced when using the combined scheme. The 
errors in vorticity are shown to reduce far below the errors encountered in an Eulerian 
calculation with a grid twice as fine in each direction. The errors in stagnation pressure 
are shown to be corrected implicitly by the correction of vorticity and velocity (through 
Crocco’s equation). 
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Two secondary flow calculations in bent pipes are performed. The secondary flow 
results from the imposition of a turbulent inlet velocity profile. Since this work deals 
with Eulerian calculations, only the tilting/stretching of the vortex lines are taken into 
account. The production of secondary flow is described, and the use of the Euler 
equations (instead of the Navier-Stokes equations) is motivated for this particular flow 
problem. The Lagrangian correction is employed as an alternative to expensive fine 
grid Eulerian calculations in order to reduce the numerical diffusion effects. As shown 
in the calculations, these effects have a strong influence on the strength and position 
of the secondary flow. Because numerical diffusion and read viscous diffusion are shown 
to result in different flow behaviors, the use of a combined Eulerian/Lagrangian scheme 
with lower numerical diffusion is required in order to assess the fluid real viscous effects 
when comparing with experiments. The limitations implied by the choice of an Eulerian 
solver are also addressed. A first attempt at including the near wall viscous effects is 
then reported as the introduction of a ‘law of the wall’ correction. Two bend geometries 
and flow cases are tested and the numerical results compared with available experiments. 

The external flow over a three-dimensioned wing is the object of the last flow case. 
The basic Eulerian and combined Eulerian/Lagrangian solutions are compared in terms 
of n um erical diffusion of the tip vortex behind the trailing edge. The numerical diffusion 
is shown to be excessive for the Eulerian solution and reduced for the combined scheme. 
Comparison with experiments is performed. The flexibility of the proposed method is 
shown by correcting the solution selectively in the tip- vortex region, and by m i n i m izing 
the CPU required for the calculation. 


13.2 Contributions 

To the author’s knowledge, the coupling of an Eulerian and a Lagrangian solution pro- 
cedures in three dimensions, enabling a correction of the Eulerian state vector based 
on Lagrangian values and aimed at reducing numerical diffusion errors, represents an 
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original contribution. The present 3-D version is based on the 2-D work of [20], but 
includes the computation of the source- terms for the vorticity not required in two di- 
mensions. The alternate upstream integration technique of the streamlines is also an 
original contribution of this work. 

The efficiency and flexibility of this new approach has been demonstrated by appli- 
cation to flow cases of different characteristics. The treated examples included steady, 
unsteady, compressible, incompressible as well as internal and external flow applica- 
tions. The combined Eulerian/Lagrangian scheme takes advantage of both the accurate 
‘elliptic’ representation of the Eulerian solution enforcing the mass requirements and 
setting the pressure field and the convection capturing capabilities of the Lagrangian 
solution. The contribution of this work is, therefore, the addition of built-in convection 
properties to a standard Eulerian solver. 

The numerical diffusion effects are quantified by comparing Eulerian solutions on 
coarse and fine grids to coarse grid Eulerian/Lagrangian solutions. Using the same 
solution comparisons, the CPU requirements are shown to be substantially lower, for a 
given accuracy, when using the combined Eulerian/ Lagrangian scheme. Also, reducing 
the numerical diffusion allows one to identify the true inviscid behavior of the flow and 
to assess the real viscosity effects when comparing with experiments. 

In the course of this work, a vorticity gradient augmentation phenomenon has been 
identified which resulted in poor vorticity representation on the Eulerian grid and desta- 
bilization of the combined scheme. In order to remediate to this problem, a Lagrangian 
pseudo- diffusion term has been added to the He lmh oltz equation, without substantially 
compromising the reduction of numerical diffusion. 

In its present form, the Lagrangian correction technique consists in a set of sub- 
routines which can be added a posteriori, not only to the present Lax-Wendroff time- 
marching technique, but to any time-marching finite- volume algorithm for the solution 
of the Euler equations on hexahedxal cells. The application of the Lagrangian algorithm 
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to tetrahedral cells, for example, would require the redefinition of the linear functions 
within the cells. 


13.3 Conclusions and recommendations for future work 

As mentioned in the introduction, the different approaches taken to efficiently resolve 
small size flow non-homogeneities result from a compromise between CPU/memory 
requirements and solution accuracy. Any attempt at improving the solution accuracy 
always ends up increasing the requirements in CPU/memory. The accuracy of the flow 
solution is linked either with grid resolution issues leading to grid refinement strategies 
or with the use of a priori known solutions and analytical flow models correcting the 
basic flow solution, or also with the use of more accurate discretization/integration 
algorithms. All of these are used separately or in conjunction with each other in order 
to cope with the poor solution accuracy encountered in finite-difference solutions of the 
Euler (or Navier- Stokes) equations in regions of high gradients and small length scales 
with respect to the background flow. All these studies were prompted by the important 
influence of these regions on the overall flow solution. 

The objective of this work was to obtain a more accurate solution on a fixed grid 
and can also be viewed as a lowering of the CPU/memory (grid size) requirements for 
a given accuracy. The comparison of the proposed approach against existing methods 
b&s been expanded upon in the introduction and will not be repeated here. 

By increasing the solution accuracy on the given grid, the proposed combined Eule- 
rian/Lagrangian scheme allows for a more accurate solution of the high gradients of the 
flow. Two of the main advantages of the proposed method are the flexibility by which 
corrections to the Eulerian solution can be performed in only chosen areas of the compu- 
tational domain, and the fact that the markers have an effect limited to the cells where 
they are located. This makes the combined Eulerian/Lagrangian scheme perform better 
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in situations where the flow non-homogeneities axe concentrated, and less suited where 
the vorticity and entropy corrections have to cover the entire domain. For example, the 
correction of the relatively small non-homogeneity represented by the unsteady convec- 
tion of a Lamb vortex in a channel (Chapter 7) leads to ~ 30% CPU increase over the 
basic Eulerian solution. In the case of steady pipe flows, since markers have been placed 
at the center of each cell, the CPU increase due to the introduction of the Lagrangian 
correction is directly linked to the size of the computational domain. The introduction 
of the Lagrangian correction for the preservation of a turbulent inlet velocity profile in a 
straight pipe (Chapter 8) results in a 30% CPU increase. The correction of the swirling 
flow of Chapter 9, in the identical computational domain, required a 33% CPU increase. 
In the case of the constant stagnation pressure flow in the 90° bend of Chapter 10, with 
twice as many nodes as the previous pipe cases, the increase in CPU was ~ 75%. Of 
course, the particular CPU increase depends also on the flow characteristics (a high 
turning of streamlines requiring more Lagrangian integration steps). 

An advantage of the combined scheme is the presence of built-in convective prop- 
erties, so that an a priori knowledge of the position or strength of the flow non- 
homogeneities is not required. 

Another conclusion to this work is that the use of the combined Eulerian/ Lagrangian 
scheme does not quite eliminate the need for extra grid resolution. In reducing the 
numerical diffusion phenomenon, the gradients of the flow become larger and, in some 
flow situations, the grid becomes too coarse to accurately support the new gradients, 
as mentioned in Section 9.2. The same grid resolution issue limits the correction of the 
wing tip vortex from only 20% chord downstream of the trailing edge (upstream of this 
distance, the vortex was defined on too few cells to attempt a Lagrangian correction 
procedure). This phenomenon is particularly acute in in viscid flow computations where 
no viscous diffusion occurs to spread the non-homogeneity on more grid cells. 

The use of the proposed Eulerian/ Lagrangian technique in combination with a grid 
adaptation should prove to be a possible solution to this problem. In addition, the 
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refinement of the grid could be linked to the markers trajectories. For example, when 
placing markers at the inlet of a bent pipe in the high vorticity regions near the wall, the 
grid could be refined along the markers trajectories as the markers move downstream 
and gather into a tight secondary vortex. The level of grid refinement could be based 
on the vorticity and entropy carried by the markers. Also, this combination technique 
would be advantageous in terms of CPU requirements for pipe flows where, instead of 
placing markers in each cell, the selective addition of markers could be realized since the 
refinement technique would prohibit many markers with largely different state vectors 
from being located in the same cell during the tightening of the secondary vortex. 

The use of a Lagrangian particle tracking solver in conjunction with the Navier- 
Stokes equations instead of the Euler equations requires the addition of the viscous 
dissipation and diffusion term in the Helmholtz equation. This is not viewed as a 
main roadblock since a similar term has already been implemented in this equation 
when using a Lagrangian pseudo- diffusion along the markers streamlines in Section 9.2. 
As mentioned by Drela in [20], the main problem here is to account for turbulence (as 
suitable models for turbulence do not exist for vortex wakes). Fortunately, the dynamics 
of vortex wakes is mostly inviscid. However, in the case of a secondary flow calculation, 
while the combined Eulerian/ Lagrangian solver gives realistic results near the entry of 
the bend (or bends of small turnings), the modeling of the migration of the viscous 
effects towards the pipe center would bring the computed results closer to experiments 
for high degrees of turning. In this case, the Lagrangian technique would not be used 
near the pipe wadis since the resolution provided by a Navier-Stokes grid ensures accurate 
solution of the viscous stresses near the walls. Still, the Lagrangian technique would 
be useful near the pipe center where the grid resolution is comparatively coarse for the 
accurate solution of the secondary vortex. 

The Lagrangian equations defined in Chapter 4 are not valid through a shock. The 
implementation of the combined Eulerian/Lagrangian scheme for shock flows is feasible 
if the markers do not cross the shock region. Therefore, the Lagrangian state vector has 
to be reinitialized (using the Eulerian solution) from one side of the shock to the other. 
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Since the correct jump conditions across the shock are ensured for any conservative Eu- 
lerian scheme, this procedure should provide the markers with correct values of entropy 
and vorticity. Also, as mentioned in [20], the exact location where the markers have to 
be stopped and reinitialized should not be critical since the usual grid refinement used 
in Eulerian solutions near the shock region will ensure a locally accurate solution. 
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Appendix A 


Mesh Generation 


The meshes used in this work are generated by solving a set of partial differential 
equations (PDE) of Poisson type using an iterative procedure, see Reference [ 75 ] for 
instance. A boundary-conforming procedure using a structured mesh is employed which 
consists of mapping the three-dimensional physical domain of Cartesian coordinates 
(z,y, z) onto a cubic computational domain of coordinates (£,?7,C)* The first step in 
defining the transformation is to specify Dirichlet and/or Neumann conditions on the 
limi ting surfaces of the physical region; those boundaries being represented by a constant 
£ or 77 or ( in the computational domain. In the elliptic partial differential method, the 
distribution of the interior grid points is then governed by the following Poisson system 

Cxx + £yy 4 " £22 — ^(£5 V 5 C) 

Vxx + Vyy + Vzz = <?(£,»?, C) 

Cxx + C yy + Ciz = R{t,V,C) t (A.l) 

where P, Q and R are source terms that can be selected to control the mesh points 
distribution. Since it is much easier to solve a system of PDE on the uniformly spaced 
grid of the computational domain, it is useful to transform system (A.l) onto the com- 
putational space. This is achieved by interchanging the roles of the dependent (£,?7,C) 
and independent (2, t/, z) variables in Eqs. (A.l). This yields an elliptic system of quasi- 
linear equations that can be written in the vector form [ 65 ] 

a n(*tt + ^) + a 22(f r T 7 r ? + ^r 7 ) + a 33 (r cc + Ar c ) + 2 (a 12 r^ + a 13 f^ + a 2 3r v(: ) = 0 , (A.2) 
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where f — (z,t/,z) is the position vector, 


dij — ^ ^ A m iA. 


mj 


m= 1 


and Ami is the cofactor of the (m, z) element in the following matrix 


M = 


The forcing functions <£, z/? and A serve to control the interior mesh points distribution, 

, J 2 P , J 2 Q . J 2 f? 

an a22 a 33 


x i 

x v 

H 

Vi 

Vv 

VC 

H 

z v 

Z C 


where 


_ d(z,3/,2) 


det^r^f^). 


Given a proper choice of the source terms P, Q and P, this transformation defines a 
one-to-one correspondence between the two spaces. The source terms Eire automatically 
evaluated in order to provide a control of the cell size and the skewness at the chosen 
domain boundaries according to the general procedure described by Sorenson [66]. 


Depending on the flow cases, two types of mesh generation strategies are used in 
this work. For internal, pipe flow computations, a 2-D version of the elliptic PDE 
solver is used to generate the cross flow grid planes, i.e. surfaces approximately normal 
to the main flow direction. These mesh planes are then arranged in the streamwise 
direction along the bend, see Figure A.l. Notice that due to the symmetry of the 
problem, only half of the pipe geometry is actually needed. The 2-D mesh plane is 
defined by a mixed 0-H topology, see Figure A.l, with the O mesh located close to the 
pipe wall, allowing for a good control of the spacing and orthogonality at this boundary. 
The mesh singularity at the center of the pipe is removed by filling this part with an 
H-type grid. The connection between the two types of grid generates two field nodes 
that are surrounded by only three cells instead of four. As discussed in Section 3.3 
this leads to a minor change in the numerical algorithm. However, this mixed grid 
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Figure A.l: Cross-sectional view of grid and distribution of cross-sections along the 
pipe. 

topology provides much more flexibility in terms of nodes distribution than a single 0- 
or H-type. Also this removes any mesh singularity from the pipe boundary and the 
symmetry plane. Notice that the two singular field nodes are allowed to move during 
the iterative solution procedure according to the number of nodes allocated for the H 
and 0 parts, and the grid control parameters set at the wall and the symmetry surface. 
In fact the position of these singular nodes is computed as an average of the position of 
the three surrounding nodes, thus enforcing identical cell volumes in this area. 

For these pipe geometries, grid control is applied on the solid wall as well as at the 
symmetry surface. This provides a good control of the discrete representation of the 
boundcCry-layer velocity profile near the wall. Grid control at the symmetry surface is 
required in order to get a good representation of the secondary vortex generated close 
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to the inside wall. This vortex is then moving along the symmetry plane due to the 
presence of its image and the wall. 

For the wing flow computation of Chapter 12, the full set of 3-D elliptic equations are 
solved to define the interior domain nodes distribution. The initial distribution of nodes 
is obtained by setting up the wing cross-sectional planes perpendicularly aroung the 
airfoil. As an option, the elliptic procedure can also be performed in two dimensions for 
each of the cross-sections, to define a smoother initial condition for the three- dimensioned 
elliptic grid solution. For this case a C-H type of mesh is used. The C part is required 
for a good resolution of the strong gradients at the leading edge and the wake flow. The 
distribution of the nodes on the wing surface is established first using multiple cubic 
spline functions along the span and the chord. This allows to cluster the nodes near 
the leading and the trailing edge as well as close to the wing tip, see Figures A. 3 and 
12.2. A set of mix ed Dirichlet/von Neumann boundary conditions are used to define the 
symmetry plane and the outer surface. Control of the spacing and the orthogonality is 
applied on the wing surface, in the wake portion (i.e. from the trailing edge to the exit) 
and also on the outer surface. The use of a symmetric profile allows for the generation 
of only one half (upper or lower) of the mesh, which is then duplicated and finally 
tilted upwards from the trailing edge to the exit to follow approximately the upward 
movement of the wake and trailing vortex, see Figure 12.3. 
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Figure A. 2: C-H mesh structure shown by different mesh surfaces. 
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Figure A. 3: Detail of wing-tip region for a streamwise cross-section located at ~ 50% 
chord with clustering in the direction perpendicular to the wall and near the tip (the 
derivatives are also prescribed using the source-terms). 
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Appendix B 

Volume and area calculations 


This appendix presents the computation of the cell volumes, the faces area and the 
volumes associated with the nodes which are required by the Lax-Wendroff algorithm. 

The volume V c of a cell is computed as the sum of the volumes of the five constituent 
tetrahedra shown on Figure B.l. 

Vc = ^5627 + F5214 + V 5 874 + V2734 + V 5274 , (B.l) 

where the indices refer to the nodes constituting each tetrahedron. The volume of a 
tetrahedron is found through a determinant constructed by the tetrahedron’s vertices. 
For example, the volume of tetrahedron 5627 is found as 

1 x5 y5 z5 
1 x6 y6 z6 

*5627 = - 

1 x2 y2 z2 
1 xl yl z7 

The area of a cell face is given by half the cross product of its diagonals. The surface 
vectors are chosen to point outwards find are written as 

51 = i(r 3 - r 6 ) x (r *7 - r 2 ) 

5 2 = -(r 8 - r*i) X (r 4 - r 5 ) 

5 3 = -(r 7 - r 4 ) x (r 3 - r 8 ) 
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(B.2) 


Si = 7^6 - n) X (^5 - ^ 2 ) 

Ss = - f s ) x (f 8 - fe) 

Se = ^3 - n) x (^2 - ^4) 

Figure B.l shows the vector surface numbering S\ to 5 6 . 

Si 



Figure B.l: Dividing of cell into five tetrahedra and surface vectors numbering definition. 

The Lax-Wendroff algorithm also uses volumes associated with nodes. A ‘node 
volume’ is defined at node n as the average of the volumes of the eight surrounding 

cells. 

= v cdl . (B.3) 

scells 

In the case of nodes lying on boundary surfaces like a wall surface or an inlet/exit 
surface, only four cells or less (comer nodes) contribute to the volume associated with 
the node. For a wall surface case, this is taken into account when applying the wall 
boundary condition as mentioned in Section 3.5. 
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Because of the symmetry condition, the volume associated with a node lying on 
a symmetry surface is found by doubling the node volume contribution from the four 
existing boundary cells. 
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Appendix C 


Stability analysis 


In order to perform a stability analysis on the Euler equations, the system is first recast 
in primitive form and computational coordinates. 


c.l Primitive form of Euler equations in computational 
coordinates 


Starting from the Euler equations in conservative form and Cartesian coordinates 


dU dF dG dH _ 
dt dx dy dz 


(C.l) 


and applying the chain rule derivation to the above equations gives 


au dF dF dF dG dG dG . dH , dH dH , n x 
dt + d^ x ar ? T?x+ ac Cl+ d(^y + dr 1 Tly+ ac Cv+ d( ^ + ~d ^ Vz+ -0, (c ' 2) 


where index x is short for derivative with respect to x, and so on for y and z. J 
is the Jacobian of the transformation from the ( x,y,z ) physical space to the (£, rj, £) 
computational space and is defined as 


J = X((y v Z( - y (Zrj ) - x v {y ( z < - y ( z ( ) + - y n Z(.). 

The metrics of the transformation are 

_ V-nH ~ y<l z n e _ x v z C ~ x < z v , x vV< ~ x (Vv 
^ - j ’ tv “ j ’ «« - j 


(C.3) 


(C.4) 
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Vx = 

ViH - y< z (. 

j 

Vy = 

x tH - X CH 
J 

x ty< - x at 

TU= j > 

(C.5) 

Cx = 

y<Zr, - y v Zi 

J 

Cy = 

x i z n ~ x v z t 

j 

^ X(Vr, - X V y< 

Cz " J 

(C.6) 


Noticing that 

■§^(yri z <: - vc z v) — -§^{y(. z <i ~ y < z t) ~ ^(yt 2 *) ~ y v z t) = ®> 

-$g( x v z <: - x c z v) ~ Jfij( x t z < ~ x c z c) ~ j%( x t Zr i ~ x v z () = o* 

■§^( x vy< - x <yv) - jftj( x zy<i ~ x <yt) ~ - x vVi) = 

the Euler equations recasted into computational coordinates are 

3U dF_ d&_ dlT 

J dt 8£ drj d( 

- f j£{yr, z < - yc z v) - -j^(y( z < ~ y < z t) ~ j£( y t Zr > ~ yr > z ^ 

=0 

- G -^{x v z C - x C z v) - ^( x i z < - X (H) - - x v z () 

=0 

■ ^ 0 ^ 

- B ^(*, , 2 /c - HVv) - fy( x ( y < - x < y i) - dC^ Vr > - = °» ( C - 7 ) 

> - - ^ y 

=0 

where F\ G\ S' are the contravariant fluxes defined as 

F' = (y v z< - y ( Zr,)F - (z„z c - x ( Zr,)G + (x v y c - x^H, 

G' = -(y^z^ - y(Z{)F + (®{Z< - Z(Z ( )G - (x^y^ - x ( y^)H, 

H' — (yfz,, — Vv z ()F — ( a: f' 2 »j — x n z i)G + — x vyc)B- 

The following relations will be used in the next section. 

r x = u{y v z ( - y(Zr,) - v{x v Z( - x^) + w(x v y ( - *<y„), 

r 2 = - u(y (Z ( - y ( zt) + v{x ( z ( - x ( z c ) - to(a:^ c - x ( y ( ), 

r 3 = u(y ( Zr, - y v z ( ) - v(x (Z t , - x v z ( ) + w(x ( y n - x„y ( ), 

whose derivatives are 

= J£( y ” z < - K**) - §jk*n* - *c^) + ^(*ny< - *c%). ( C - 8 ) 
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dr 2 du dv . dw . . . 

= ~^(y( z < - y< z t) + ^( x ( z < - *< z c) - - *cvc), (c. 9 ) 

dr 3 du dv . dw 

~ iy( z v — yv z () ~ ~q£' x ( Zt > ~ x v z i) + ~Q^\ x iyv — *•?%)• (c.io) 


C.1.1 Primitive form: incompressible flow 


In the case of an incompressible flow, the state vector U and the contravariant fluxes 
F',G',H' are 


U = 


G' = 


/ 

r 2 u - p*(y i z < - y ( z ( ) 
r 2 v + p*(x£Z( - x ( Z() 

\ r 2 w -p*(x ( y ( - z ( y ( ) ) 

The continuity equation is 

r dp* 


(p') 



u 

, f = 

r i u + p*(y v z( -yaz n ) 

V 


r lV -P*{x v z c - x ( z v ) 

{ w ) 


r\w - p*(i n y c - j 

c\r 2 

\ 

1 ❖» 


, H' = 


r 3 « + P*(»^ ~ VrjZi) 

- v*{ x i z n - x v z t) 

V r 3 W ~ P*( X ty r, - X r,V() ) 


T ~r 2 (9ri dv 2 dr z \ 

dt ° V dt dv + dc ) ’ 


(C.ll) 


where the RHS is found from Equations (C.8) to (C.IO). The first component of the 
momentum gives 


_ du 

J ~di = 


du du du ( dr i dr 2 dr 3 \ 

r ' H " T, d~r, - rs K - “ V af + a? + a<j 

(Vn z C - yt 2 *)^ + (V(H - - (*(** - 

( d d d \ 

o^(yv z < - y< z v ) - fyiy^ - yen) + y - vvh)) • (c.12) 


=0 
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Similarly, one finds for the two other components of the momentum equation 



dv dv dv 

Tl l^ ~ T2 dv ~ r3 ~dC ~ V 


(dj\ f^2 dr 3 \ 

\ d£ dr) d( ) 


+ 


( x v z < - x d z n)^ ~ ( X € Z C 


\ dp* < \ dp* 

X <. Z i)~Q^ + ( X t Zr > ~ Xr > Z () ' 


(C. 13 ) 



OW OW OW ( OTi OT 2 or 3 \ 

ri ~dt ~ T2 !hj ~ T3 !k ~ w V ae + + ac ) 

(x v y< - x <vv)-j£ + ( x tvc - ~ ( x t yr > ~ Xr>y ^~§£’ ( C-14 ) 


Collecting Equations (C.ll) to (C. 14 ) gives the Euler equations modified by the artificial 
compressiblity concept and written in primitive form and computational coordinates 


t 


du T -i f a du . du . 

Ht +J i'W '~Sv 1 


sc.) 


where 


U = 


u 

V 


\ w ) 


1 o c l{vv z c - y< z v ) ~ c l( x v z < - x a z n ) c 2 (*r,y< - x <y n) N 

(y^c - yc 2 *?) u ivri z < - y<: z v) + r i ~ u ( x v z c - x c z v) u ( x vy< - x cVv) 

- X^Zr,) v(yr,Z( - y C Zr,) -V ( Xr ,Z ( - X^) + Ty v(x v y ( - * ( y„) 

y (x„y< - x c yr,) - y c z v ) -w(xr,z ( - x^) w{x v y < - x ( y v ) + t x J 


■ o -cliviH - y< z t) c K x ( z c - x ch) - c K x (y< - x < y ( ) ^ 

-(y<z< -y<z t ) -u{ytZ( - y<z t ) + r 2 t i{x t z c - x^) -u{x iy< -x <y*) 

(x i z c - x ( z f ) -t>(y f z c - y c z f ) v{x^ - x<z f ) + r 2 -v{x^y^ - x^) 

\~( x cyc ~ *( 2 It) -w(y{ 2 ( - y^zO w{x^ - x ( z ( ) -w(x ( y ( - x <yf) + r 2 ) 


1 0 cliVt 2 * ~ VvH) — “ x v z i) c l( x tVv ~ x vVt) ^ 

(y«Zr T~yn z t) u (yt z ri ~ Vv z t) + r 3 -u(x ( Zr,-x r , z () «(*<» ^~ x vV() 

-(x^Zr, - XrjZi) v(y ( Zr, - JfaZf) -«(*{*, - ® n Z() + r 3 

- air,2/«) -»(***!-*•»**) U, ( a: fy») - *r,y«) + T v 
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C.1.2 Primitive form: compressible flow 


In the case of a compressible flow the state vector U and the contravariant fluxes 
F\ G', H' are 


pr lU + piyrjZt - y c z v ) 

, F' = p r lV _ p(x r ,z < - x ( Zr,) , G' = 
pr\W + p{x v y c - x ( y v ) 
l ri{pe 0 +p) i 


pr 2 u - p(y i z < - y c z ( ) 
pr 2 v + p{x{Z( - X(Z() , 

pr 2 w - - x^) 

r2{peo + p) , 


pr z u + p{y^z n -yr,Z() 
H' = pr 3 V - pix^Zr, - Xr,z{) 
pr 3 w+p{x ( y v - Xr,y ( ) 
i r 3 (pe 0 + p) j 


The continuity equation is 


jdp d d d 

J Tt = 


dp dp dp 

~ ~ r 'T( ~ r ’^ ‘ ~ p 


dr\ dr 2 dr 3 

lit + Ihi + ~d( 


The first component of momentum is 


du ( d . ( 

Jp Tt = 3 U (pu) - u ] 


(C.15) 


(C.16) 


Replacing the expression for d(pu)/dt and using Equation (C.15) for dpfdt , the first 
component of momentum becomes 

du du du du dp, dp, x dp, x 

Jp Tt = 

(0-17) 
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S imil arly, one finds for the two other components of the momentum equation 


dv dv dv dv dp, * dp, . dp. . 

(C.18) 

dw dw dw dw dp. . , dp. . dp. , 

(C.19) 

By t akin g the time-derivative of the equation of state we get 

1 5p , (d , 1. , , dp du dv dw\ 

J T^Tt = 3 - 2 ( “ + * + ” * m - Tt - ~ '"»? ) ' (C20) 

Replacing the expressions for d(pE)/dt and using Equations (C.15) and (C.17), (C.18) 
and (C.19) we get 

J?P = _ ri ^ _ r jL - r 3 °P - 7 p ^ (C.21) 

dt d£ 2 dr} 3 dC \dt dp d() 

Thus, the Euler equations for a compressible flow written in primitive form and in 

computational coordinates are 




dU 


dt 


dV ■ 


dU v 


dU T 


p + J- 1 Ia 1 ^ + a 2 ^ + a 3 ^i> 


dt 


dp 


dC 


)=o, 


where 


u p = 


U 

V 

w 


\ p / 


A , = 


fi p{y v zc - y< z v ) -pi x v z c - p( x M - x <y*i) 0 

Or! 0 0 i(yr,z c - y^Zr,) 

0 0 r : 0 -^{x n z ( - 

oo o n ${* v yc - x <y * ?) 

o ip{y, n z( - yc z v) -ipi x v z c - x c z v) ip( x r,y< - x cVv) r i 
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C.2 Stability 


The discretization following Ni’s method gives 

Up +1 = Up + (AifitUrfSt + A^nc^Sr, + A 3 Hr,(J-i6c) U£ + (C.22) 

^At 2 (Ai^UrjS^ + A 2 fi^(6 v + JhMrj^c) 2 
where the operators 6 and /z axe defined by 


6(Uijk — 

Ui+%jk - Ui-Ljk' ^Uijk - U lj+ i k U { -_i k , 

(C.23) 

6(Uijk = 

u ijk + i - u ijk-\* 

(C.24) 

H(Uijk = 

2^i+%jk + 

(C.25) 

H<Uijk = 

|(>Wt + 

(C.26) 
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Considering a wave-like solution U£ = the amplification matrix G 

is 


G{6i, 02, 0 3 ) = I - 2 i COS Y COS y COS yAi23 + 2A f 2 A 2 23 , 


(C.27) 


where 


A A . 01 #2 i A * 03 4 . 03 01 02 

^123 = A 1 Sin - COS y COS y + A 2 sin — COS y COS y + Sin y COS y COS y . 


Defining 


(C.28) 


D — At ($i Ai + s 2 A 2 + S3A3), 


(C.29) 


// • 01 02 03 v 2 / • 02 03x 2 . / 02 . 03 \ 2 ( n on\ 

5 = v ( sm y cos y cos y ) + ( cos y sm y cos y ) 2 + (cos — cos y sm — ) 2 , (C.30) 


sin 4- cos cos % 


cos ^ sin % cos ^ 


$1 — 


5 2 = 


^3 = 


cos ^ cos ^ sin ^ 


(C.31) 


with the property sj + sj + s| — 1, the amplification matrix G can be written as 


(?(0i,02,03) = J — 2 t s cos y cos ^ cos y J 9 - 2 s 2 D 2 . 

Z Z / 


(C.32) 


If Ap is an eigenvalue of D, the corresponding eigenvector is also an eigenvector of G 
with an associated eigenvalue 


a 1 O ' & 1 03 \ o 2 \ 2 

Ac = 1 — 2z s cos — cos COS — AD — 2s Af) . 

z z z 


(C.33) 


If Ap is reed and |Aj>| < 1, then 
|Ag | 2 = 1 — 4s 2 A£j i\ — s 2 X 2 d - cos 


01 2 02 2 03 2 \ 

COS — cos I 

2 2 2 y 

t) 


01 2 02 2 


< 1 — 4s 2 Ap 1 — s 2 — cos ^r- cos ~~ cos 


/ id 2 a 2 /j 2 a 2 /i 2 /j 2 a 2 a 2 n 

i a 2 \2 ( • “l • "2 P 3 . ^1 . . &3 . . &1 “2 . ^3 

- 1 — 45 An sin — sm — cos — + cos — sm — sm — +sm — cos — sm — 

D \ 2 2 2 2 2 2 2 2 2 


< 1. 
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Thus the condition X D < 1 is sufficient (but may be not necessary) for stability. 


C.2.1 Stability: incompressible flow 


Let 


r = riSi + r 2 s 2 + r 3 s 3 , (C.34) 

a = {y v Z{ - y^Zr,^ - (ytz ( - y ( z i )s 2 + (y^Zr, - y^z^ss, (C.35) 

b = ~(x ri Z <: - X(Zr,)si + [X(Z( - X(Z{)s 2 - ( x t z V - x v Z t) S 3 > (C.36) 

d = {x v y( - x c yt)si - (x^ - X(y()s 2 + {x^y v - x v y c )s 3 . (C.37) 

Using Equation (C.29) and the previous definitions for r,a, 6 and d, the matrix D is 
written as 


/ 


D = 


cld 


0 c 2 a a c 2 b 

a ua + r ub ud 

b va vb + r vd 

d wa wb wd + r J 


The eigenvalues of D are 

^ (r, r, r + r - ^ + c2(a 2 + 6 2 + d 2 )) . (C.38) 

A conservative estimate for the maximum eigenvalue is 

A Dm „ < Y (>/ r l +r 2 + r 3 + >/ r l + ** 2 + r 3 + C a(« 2 + ^ + ^)) ’ ( C ‘ 39 ) 

where the value of r has been maximized by yjr\ + t\ + r| since sf + s 2 + = 1, and 

the values of a 2 ,b 2 ,<P have been maximized by a 2 , b 2 , <? with 


a 2 = (y v z ( - VcZv) 2 + {V( z < ~ V< z t) 2 + ~ ' (C.40) 

b 2 = {x v Z C - X ( Zr ,) 2 + {x C Z< - X(Zt ) 2 + {xtZr, - X^) 2 , (C.41) 
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d = ( x vV< - x cvc ) 2 + ( X (V< - x CVt ) 2 + ( x CVv - Xr,yt:) 2 . (C. 42 ) 

Thus a sufficient criterion for stability is 

At < — L . (C.43) 

\J r| + r\ + t\ + \j r\ -f t\ + r§ + c 2 (a 2 + b 2 +<?) 


C.2.2 Stability: compressible flow 


From Equation (C. 29 ) and using the definitions for r, a, b and d, the matrix D is 


r pa 
0 r 


D = 


0 0 
0 0 


pb pd 0 

0 0 a / p 

r 0 bf p 

0 r d/ p 


y 0 7 pa 7 pb 7 pd r J 


The eigenvalues of D are 



r, r, r + c\J a 1 + b 2 + d 2 , 


t — cVa? + b 2 + d 2 j , 


(C. 44 ) 


where c is the speed of sound c = y/rpfp* A conservative estimate for the maximum 
eigenvalue is 

^ (\Ai 2 + r ! + r 3 + cv/a 2 + fc 2 + d 2 ) , (C. 45 ) 

where the value of r has been maximized by yjr\ + r| + r| since s 2 + + s| = 1. Thus 

a sufficient criterion for stability is 

At < - . (C. 46 ) 

yjr \ + t\ + r| + c\/a 2 + 6 2 + d 2 
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Appendix D 


Brute force location of markers 


This type of marker location in the whole domain is used for the initial location of the 
markers or as a last resort search when the marker has moved over a distance larger 
than one cell from its previous position and the usual search in the neighbouring cells 
has failed. First, nearly all cells are disqualified by eliminating the cells whose ranges 
of nodes coordinates do not encompass the marker coordinates. Each remaining cell is 
then divided into 6 tetrahedra T1245, T 456 8, T 24 56> T 2 346 5 TWs, T 36 78 (the indices refer to 
the node numbering defined in Appendix B). A search is performed in each tetrahedron 
by expressing the marker position vector p defined on Figure D.l in the local coordinate 
system (ei,e2,e3) as p = ae\ + /3?2 + 7^3 where a, 7 are found by the relations 

(e 2 x €3) • p (ei x £3) • p (eixe 2 )-p 

(«i x^)'^’ (e\ xejJ-eV (eiXe 2 )-e 3 

The necessary conditions for the marker to be located in the tetrahedron are 

«>0, /3 > 0, 7 > 0, a + 0 + 7<1. (D.2) 



Figure D.l: Loccil coordinates in tetrahedron. 
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Appendix E 


Newton-Raphson procedure 


E.l Marker location in cell 


The (£,r),Q coordinates of a marker located at ( x,y,z ) in a given cell are found by 
solving the implicit system 

8 

r = £^(£,77,0^ (E.l) 

fc=l 

by a few Newton iterations. Or, in shorthand notation, 

9 2(£,*?,C) 

_ 77, C) 

Starting from an initial guess Qo and using a Newton-Raphson iteration, the location 
of the marker in (f, 77, £) is found as 


= <?(«) = r - £ N k (a) f k = 0 . 


(E.2) 


fc=i 


ff(a) = 9(^0) + 


dg 


da 


Aa. 


(E. 3 ) 


J C*o 


The Newton-Raphson procedure consists then in solving the following equations itera- 
tively until convergence is obtained, that is when Aa — ► 0 . 


r -1 


Aa = — 


dg 


da 


g{a 0 ), a 0 = a 0 + Aa 


(E. 4 ) 


-* cr 0 
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The first guess for a 0 is set to 0 so that the search begins from the cell center. In practice 
only three iterations suffice to determine the location of the marker with acceptable 
accuracy. The Jacobian matrix is obtained from Equation (E.2) as 

% = - 1 < E - 5 > 


fc=i 


da 


Using Equations (5.1) for the definition of the tri-linear functions Ni to N% gives 


|| = | [(1 - 77)((1 - C)(n - ^ 2 ) + (1 + 0(^5 - ^e)) 

+ (1 + r/)((l - C)(r 4 - r 3 ) + (1 + C)(r*8 - r 7 ))], 

f| = J [(1 - 0((1 - C)(ri - u) + (1 + C)(r*s - r s )) 

+(1 + 0((! ~ C)(r *2 - r*3) + (1 + C)(r*6 - r 7 ))], 

|| = 5 [(1 - 4)((1 - 7 ?)( f i ~ f s) + (1 + V){U ~ r 8 )) 

+(1 + 0 ((! “ 7 ?)(^2 - r 6 ) + (1 + T 7 )(r 3 - r 7 ))]. 


E.2 Metrics derivatives at marker location 


The terms (da/ dr) required for the evaluation of the derivative of the function / at the 
marker location in Equation (5.8) are obtained by deriving Equation (E.2) as 


Hence 


da 

dr 



d£ 

dr 



da 

= ~dS' 


4y 






[" dr' 

1 \dg' 

Vv 

dz 

~ [ds. 

[da. 


(E. 6 ) 


(E.7) 


[ Cx Cy Cz J 

where the indices a 16 shorthand for d/d£,d/dr],d/dC. Or explicitly 


, g2 v gZc - 9^c9^v f _ sM 3 c - lk?h c - gM 2 c - g 1 cg 2 >> 

4* = — 1 — -j > 4y J > - 7 » 

g2 e gZ c - g2 ( g3i _ gl f g3 c - gl ( g3 { _ _ 9^j - g^g gj 

Vx ~ J 1 Vy — J ’ 'i z J 5 

, g2fg3 n — y2 n </3( , _ glfg3 n ~ 9bgZj f _ glfg 2 n ~ gM 2 f 

Cx = — J > Cy > « - 7 ’ 
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where 


J = -gl^g^ - g2 c g3 T) ) + glr,(g2 ig 3 ( - g 2 <g Z ( ) - gl ( (g2 ( g3 v - g2 T ,g3 i ). (E.8) 
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Appendix F 


Change in Circulation Due to Diffusion 


The loss in circulation due to diffusion can be expressed by writing the convective 
change in circulation around a closed curve as 

157 = dr) = / dr) = f^' dr = f~dt' dr ' (F,1) 

±§v-v=0 

Using the momentum equation to substitute for the convective change in circulation 
is 

£L = j(-^).dr+ j(u^ 2 v)-dr. (F.2) 

By using Stokes theorem, the first integral is expressed as 

^(_^).d r = ^(-V(i)x Vp)-d5, (F-3) 

which gives a zero contribution for incompressible flows. Transforming the second inte- 
gral of Equation (F.2) with 


V 2 r = V x (V x v) + V(V^_v) = -V x w, (F.4) 

=0 

the convective change in circulation due to viscous effects in an incompressible flow is 
written as 


Dr 

Dt 



(F.5) 
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